Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add a GDAL VRT dataset descriptor for the SRTM datasets #16

Closed
Ryanf55 opened this issue Jan 10, 2024 · 6 comments
Closed

Add a GDAL VRT dataset descriptor for the SRTM datasets #16

Ryanf55 opened this issue Jan 10, 2024 · 6 comments

Comments

@Ryanf55
Copy link

Ryanf55 commented Jan 10, 2024

The SRTM dataset hosted by the ArduPilot terrain server is useful for more than just ground stations. Companion computers can also make use of the terrain data for path planning algorithms.

As part of some ongoing work to enable SmartRTL for planes, which needs to run on a companion computer, we're using GDAL as a library to load terrain data. GDAL has support for many file formats, including the terrain hosted on ArduPilot's server.

What I would like to request is to add a VRT dataset to support GDAL-based virtual file access of the terrain data.

See ethz-asl/grid_map_geo#35 for some more info and examples on what the VRT dataset looks like.

Assuming the process I have running completes, would it be possible to upload the resultant .vrt file to the AP terrain database? It makes more sense for this to be hosted by AP, since

  1. It takes a long time to compute the VRT database
  2. The file is plaintext and not very large; it's just references to the existing datasets
  3. It is useful for anyone using GDAL who wants access to AP's terrain data.

For example, I would to add a new file to this directory:
https://terrain.ardupilot.org/SRTM1/

It would be called srtm1.vrt, and it would contain the 23997 hosted files.

I took a look at what QGC implemented, and unfortunately, it is not abstracted in a way that can be used elsewhere due to its heavy dependence on the UI framework QT.

@Ryanf55
Copy link
Author

Ryanf55 commented Jan 10, 2024

The following script generates a list of SRTM1 files hosted by ArduPilot.

#!/usr/bin/env python3

# This file runs a query on ArduPilot's terrain server to list out all the SRTM1 filenames.
# It then writes them in a file to be used by GDAL's gdalbuildvrt.
# https://gdal.org/programs/gdalbuildvrt.html

# Usage:
# Get a list of all SRTM1 files hosted by ArduPilot
# $ python3 fetch_ap_listings.py
# $ head -n 3 SRTM1_list.txt
# >>> /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/N00E006.hgt.zip/N00E006.hgt
# >>> /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/N00E009.hgt.zip/N00E009.hgt
# >>> /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/N00E010.hgt.zip/N00E010.hgt
# >>> /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/N00E011.hgt.zip/N00E011.hgt

# Build the VRT. This takes forever (hours) because there are 23,997 files.
# $ gdalbuildvrt -input_file_list SRTM1_list.txt ap_srtm1.vrt
# >>> 0...10...20...30..

import requests
from bs4 import BeautifulSoup

urls = {"SRTM1": "https://terrain.ardupilot.org/SRTM1/"}


def main():
    html = requests.get(urls["SRTM1"]).text
    soup = BeautifulSoup(html)

    with open("SRTM1_list.txt", mode="w") as file:
        for a in soup.find_all("a", href=True):
            ref = a["href"]
            if ref.endswith(".hgt.zip"):
                ref_name = ref[: -(len(".zip"))]
                file.write(f"/vsizip/vsicurl/{urls['SRTM1']}{ref}/{ref_name}\n")


if __name__ == "__main__":
    main()

The resultant VRT file can be generated (slowly) with the following:

gdalbuildvrt -input_file_list SRTM1_list.txt ap_srtm1.vrt

Once complete, the file attached is the result. It's 12MB.

ap_srtm1.zip

You can see the value of this; with an off-the-shelf tool, I can now get the altitude of any location within the dataset, and not have to deal with the way AP named the files, and it gracefully handles void data or missing tiles.

$ gdallocationinfo -wgs84 ap_srtm1.vrt  149.159350 -35.363272  
Report:
  Location: (1184974P,429708L)
  Band 1:
    <LocationInfo><File>/vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/S36E149.hgt.zip/S36E149.hgt</File></LocationInfo>
    Value: 601

@stephendade
Copy link
Collaborator

Hi @Ryanf55,

Assuming the process I have running completes, would it be possible to upload the resultant .vrt file to the AP terrain database? It makes more sense for this to be hosted by AP

Happy to add the file when you've got it ready.

@Ryanf55
Copy link
Author

Ryanf55 commented Jan 15, 2024

Hi @Ryanf55,

Assuming the process I have running completes, would it be possible to upload the resultant .vrt file to the AP terrain database? It makes more sense for this to be hosted by AP

Happy to add the file when you've got it ready.

Awesome. See the post above. It's the attached zip file ap_srtm1.zip and thoroughly tested locally. Please upload it in the current form (zipped) otherwise it's slow to download.

@stephendade
Copy link
Collaborator

Done. File's up at https://terrain.ardupilot.org/SRTM1/ap_srtm1.zip

@Ryanf55
Copy link
Author

Ryanf55 commented Jan 15, 2024

Done. File's up at https://terrain.ardupilot.org/SRTM1/ap_srtm1.zip

Awesome. It works. If you are on Ubuntu, try installing gdal-bin.

$ gdalinfo /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/ap_srtm1.zip/ap_srtm1.vrt
...
       /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/S84W179.hgt.zip/S84W179.hgt
       /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/S84W180.hgt.zip/S84W180.hgt
Size is 1296001, 604801
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000138888888898,84.000138888888884)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (-180.0001389,  84.0001389) (180d 0' 0.50"W, 84d 0' 0.50"N)
Lower Left  (-180.0001389, -84.0001389) (180d 0' 0.50"W, 84d 0' 0.50"S)
Upper Right ( 180.0001389,  84.0001389) (180d 0' 0.50"E, 84d 0' 0.50"N)
Lower Right ( 180.0001389, -84.0001389) (180d 0' 0.50"E, 84d 0' 0.50"S)
Center      (   0.0000000,  -0.0000000) (  0d 0' 0.00"E,  0d 0' 0.00"S)
Band 1 Block=128x128 Type=Int16, ColorInterp=Undefined
  NoData Value=-32768

I'll file another ticket when the other files are ready. I'll likely need direct access to the server to run the processing.

@Ryanf55 Ryanf55 closed this as completed Jan 15, 2024
@Ryanf55
Copy link
Author

Ryanf55 commented Nov 6, 2024

FYI, to those trying to use it from GDAL

gdalinfo /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/ap_srtm1.zip/ap_srtm1.vrt
gdallocationinfo -wgs84  /vsizip/vsicurl/https://terrain.ardupilot.org/SRTM1/ap_srtm1.zip/ap_srtm1.vrt 149.165224 -35.363056

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants