diff --git a/sardem/conversions.py b/sardem/conversions.py index ecccf4f..1c7920d 100644 --- a/sardem/conversions.py +++ b/sardem/conversions.py @@ -3,17 +3,13 @@ import shutil import subprocess -import requests - from . import utils logger = logging.getLogger("sardem") -URL_EGM08 = "http://download.osgeo.org/proj/vdatum/egm08_25/egm08_25.gtx" -URL_EGM96 = "http://download.osgeo.org/proj/vdatum/egm96_15/egm96_15.gtx" -EGM_URLS = { - "egm96": URL_EGM96, - "egm08": URL_EGM08, +EPSG_CODES = { + "egm96": "5773", # https://epsg.io/5773 + "egm08": "3855", # https://epsg.io/3855 } EGM_FILES = { "egm96": os.path.join(utils.get_cache_dir(), "egm96_15.gtx"), @@ -22,23 +18,25 @@ def egm_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True, geoid="egm96"): - """Convert a DEM with a EGM96 vertical datum to WGS84 heights above ellipsoid""" + """Convert a DEM with a EGM96/2008 vertical datum to WGS84 heights above ellipsoid""" if output is None: ext = os.path.splitext(filename)[1] output = filename.replace(ext, ".wgs84" + ext) - egm_file = EGM_FILES[geoid] - if not os.path.exists(egm_file): - download_egm_grid(geoid=geoid) + code = EPSG_CODES[geoid] + # Source srs: WGS84 ellipsoidal horizontal, EGM geoid vertical + s_srs = '"epsg:4326+{}"'.format(code) + # Target srs: WGS84 horizontal/vertical + t_srs = '"epsg:4326"' # Note: https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-tr # If not specified, gdalwarp will generate an output raster with xsize=ysize # We want it to match the input file xsize, ysize = _get_size(filename) cmd = ( - 'gdalwarp {overwrite} -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids={egm_file}" ' - '-t_srs "+proj=longlat +datum=WGS84 +no_defs" -of ENVI -ts {xsize} {ysize} ' + 'gdalwarp {overwrite} -s_srs {s_srs} -t_srs {t_srs}' + ' -of ENVI -ts {xsize} {ysize} ' ' -multi -wo NUM_THREADS=4 -wm 4000 {inp} {out}' ) cmd = cmd.format( @@ -47,7 +45,8 @@ def egm_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True, geoid="eg overwrite="-overwrite" if overwrite else "", xsize=xsize, ysize=ysize, - egm_file=egm_file, + s_srs=s_srs, + t_srs=t_srs, ) logger.info(cmd) subprocess.run(cmd, check=True, shell=True) @@ -103,32 +102,3 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_bounds=False): if not using_gdal_bounds: # Now shift back to the .rsc is pointing to the middle of the pixel utils.shift_rsc_file(rsc_filename, to_gdal=False) - - -def download_egm_grid(geoid="egm96"): - if geoid == "egm96": - url = URL_EGM96 - elif geoid in ("egm08", "egm2008"): - url = URL_EGM08 - else: - raise ValueError("Unknown geoid: {}".format(geoid)) - - egm_file = EGM_FILES[geoid] - if os.path.exists(egm_file): - logger.info("{} already exists, skipping.".format(egm_file)) - return - - 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 diff --git a/sardem/cop_dem.py b/sardem/cop_dem.py index cb333d5..89db3c1 100644 --- a/sardem/cop_dem.py +++ b/sardem/cop_dem.py @@ -1,5 +1,4 @@ import logging -import os from copy import deepcopy import requests @@ -9,7 +8,6 @@ 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) @@ -29,7 +27,6 @@ def download_and_stitch( from osgeo import gdal gdal.UseExceptions() - geoid = "egm08" # TODO: does downloading make it run any faster? # if download_vrt: # cache_dir = utils.get_cache_dir() @@ -38,16 +35,13 @@ def download_and_stitch( # make_cop_vrt(vrt_filename) # else: vrt_filename = "/vsicurl/https://raw.githubusercontent.com/scottstanie/sardem/master/sardem/data/cop_global.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" - + code = conversions.EPSG_CODES["egm08"] + s_srs = '"epsg:4326+{}"'.format(code) + t_srs = '"epsg:4326"' xres = DEFAULT_RES / xrate yres = DEFAULT_RES / yrate @@ -68,6 +62,7 @@ def download_and_stitch( # Used the __RETURN_OPTION_LIST__ to get the list of options for debugging logger.info("Creating {}".format(output_name)) + logger.info("Fetching remote tiles...") cmd = _gdal_cmd_from_options(vrt_filename, output_name, option_dict) logger.info("Running GDAL command:") option_dict["callback"] = gdal.TermProgress