Merge pull request #5 from scottstanie/vrt-cop-dem

Vrt cop dem

sardem/cli.py

@@ -111,10 +111,10 @@ def get_cli_args():
         help="Source of SRTM data (default %(default)s). See README for more.",
     )
     parser.add_argument(
-        "--keep-egm96",
+        "--keep-egm",
         action="store_true",
         help=(
-            "Keep the DEM heights as geoid heights above EGM96. "
+            "Keep the DEM heights as geoid heights above EGM96 or EGM2008. "
             "Default is to convert to WGS84 for InSAR processing."
         ),
     )
@@ -170,6 +170,6 @@ def cli():
         args.data_source,
         args.xrate,
         args.yrate,
-        args.keep_egm96,
+        args.keep_egm,
         output,
     )

sardem/conversions.py

@@ -71,7 +71,7 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96"):
 
     Overwrites file, requires GDAL to be installed
     """
-    if not _gdal_installed_correctly():
+    if not utils._gdal_installed_correctly():
         logger.error("GDAL required to convert DEM to WGS84")
         return
 
@@ -111,19 +111,25 @@ def download_egm_grid(geoid="egm96"):
         logger.info("{} already exists, skipping.".format(egm_file))
         return
 
-    logger.info("Downloading from {} to {}".format(url, egm_file))
+    size = _get_file_size_mb(url)
+    logger.info("Performing 1-time download {} ({:.0f} MB file), saving to {}".format(url, size, egm_file))
     with open(egm_file, "wb") as f:
         resp = requests.get(url)
         f.write(resp.content)
 
 
+def _get_file_size_mb(url):
+    # https://stackoverflow.com/a/44299915/4174466
+    return int(requests.get(url, stream=True).headers['Content-length']) / 1e6
+
+
 def shift_dem_rsc(rsc_filename, outname=None, to_gdal=True):
     """Shift the top-left of a .rsc file by half pixel
 
     See here for geotransform info
     https://gdal.org/user/raster_data_model.html#affine-geotransform
     GDAL standard is to reference a raster by its top left edges,
-    while often the .rsc for SAR focusing is using the middle of a pixel.
+    while the .rsc for SAR focusing might use the middle of a pixel.
     `to_gdal`=True means it moves the X_FIRST, Y_FIRST up and left half a pixel.
     `to_gdal`=False does the reverse, back to the middle of the top left pixel
     """
@@ -152,13 +158,3 @@ def shift_dem_rsc(rsc_filename, outname=None, to_gdal=True):
         f.write(loading.format_dem_rsc(rsc_dict))
 
 
-def _gdal_installed_correctly():
-    cmd = "gdalinfo --help-general"
-    # cmd = "gdalinfo -h"
-    try:
-        subprocess.check_output(cmd, shell=True)
-        return True
-    except subprocess.CalledProcessError:
-        logger.error("GDAL command failed to run.", exc_info=True)
-        logger.error("Check GDAL installation.")
-        return False

sardem/cop_dem.py

@@ -0,0 +1,118 @@
+# Based on https://github.com/MarcYin/Copernicus_GLO_30_DEM_VRT
+import logging
+import os
+from copy import deepcopy
+
+import requests
+
+from sardem import conversions, utils
+
+TILE_LIST_URL = "https://copernicus-dem-30m.s3.amazonaws.com/tileList.txt"
+URL_TEMPLATE = "https://copernicus-dem-30m.s3.amazonaws.com/{t}/{t}.tif"
+DEFAULT_RES = 1 / 3600
+
+logger = logging.getLogger("sardem")
+utils.set_logger_handler(logger)
+
+
+def download_and_stitch(
+    output_name, bounds, keep_egm=False, xrate=1, yrate=1, download_vrt=False
+):
+    """Download the COP DEM from AWS.
+
+    Data comes as heights above EGM2008 ellipsoid, so a conversion
+    step is necessary for WGS84 heights for InsAR.
+
+    References:
+        https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198
+        https://copernicus-dem-30m.s3.amazonaws.com/readme.html
+    """
+    from osgeo import gdal
+
+    gdal.UseExceptions()
+    geoid = "egm08"
+    if download_vrt:
+        cache_dir = utils.get_cache_dir()
+        vrt_filename = os.path.join(cache_dir, "copernicus_GLO_30_dem.vrt")
+        if not os.path.exists(vrt_filename):
+            make_cop_vrt(vrt_filename)
+    else:
+        vrt_filename = "/vsicurl/https://raw.githubusercontent.com/scottstanie/sardem/vrt-cop-dem/sardem/data/copernicus_GLO_30_dem.vrt"  # noqa
+    egm_file = conversions.EGM_FILES[geoid]
+    if not os.path.exists(egm_file):
+        conversions.download_egm_grid(geoid=geoid)
+
+    if keep_egm:
+        t_srs = s_srs = None
+    else:
+        s_srs = "+proj=longlat +datum=WGS84 +no_defs +geoidgrids={}".format(egm_file)
+        t_srs = "+proj=longlat +datum=WGS84 +no_defs"
+
+    xres = DEFAULT_RES / xrate
+    yres = DEFAULT_RES / yrate
+
+    # access_mode = "overwrite" if overwrite else None
+    option_list = gdal.WarpOptions(
+        format="ENVI",
+        outputBounds=bounds,
+        dstSRS=t_srs,
+        srcSRS=s_srs,
+        xRes=xres,
+        yRes=yres,
+        outputType=gdal.GDT_Int16,
+        resampleAlg="bilinear",
+        multithread=True,
+        warpMemoryLimit=5000,
+        warpOptions=["NUM_THREADS=4"],
+        options='__RETURN_OPTION_LIST__'  # To see what the list of cli options are
+    )
+
+    logger.info("Creating {}".format(output_name))
+    cmd = _gdal_cmd_from_options(vrt_filename, output_name, option_list)
+    logger.info("Running GDAL command:")
+    logger.info(cmd)
+    gdal.Warp(output_name, vrt_filename, options=option_list)
+    return
+
+
+def _gdal_cmd_from_options(src, dst, options):
+    opts = deepcopy(options)
+    for idx, o in enumerate(opts):
+        # Wrap the srs option in quotes
+        if o.endswith("srs"):
+            opts[idx + 1] = '"{}"'.format(opts[idx + 1])
+    return "gdalwarp {} {} {}".format(src, dst, ' '.join(opts))
+
+
+def make_cop_vrt(outname="copernicus_GLO_30_dem.vrt"):
+    """Build a VRT from the Copernicus GLO 30m DEM COG dataset
+
+    Note: this is a large VRT file, ~15MB, so it can take >30 minutes to build.
+    """
+    from osgeo import gdal
+
+    gdal.UseExceptions()
+
+    tile_list = get_tile_list()
+    url_list = _make_url_list(tile_list)
+    vrt_options = gdal.BuildVRTOptions(
+        resampleAlg=gdal.GRIORA_NearestNeighbour,
+        outputBounds=[-180, -90, 180, 90],
+        resolution="highest",
+        outputSRS="EPSG:4326+3855"
+    )
+    logger.info("Building VRT {}".format(outname))
+    vrt_file = gdal.BuildVRT(outname, url_list, options=vrt_options)
+    vrt_file.FlushCache()
+    vrt_file = None
+
+
+def get_tile_list():
+    """Get the list of tiles from the Copernicus DEM 30m tile list"""
+    logger.info("Getting list of COP tiles from ", TILE_LIST_URL)
+    r = requests.get(TILE_LIST_URL)
+    return r.text.splitlines()
+
+
+def _make_url_list(tile_list):
+    return [("/vsicurl/" + URL_TEMPLATE.format(t=tile)) for tile in tile_list]