From 151eb20b50abb9591ce676cdd634689707d4f399 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 10:44:52 -0400 Subject: [PATCH 01/11] enhancement: make integer bbox align with tiles without a shift, user needs to manually add a 0.000138 difference to all bbox numbers --- sardem/cli.py | 47 +++++++++++---------- sardem/constants.py | 2 + sardem/cop_dem.py | 8 ++-- sardem/dem.py | 79 +++++++++++++++++++++--------------- sardem/tests/test_dem.py | 16 +++++--- sardem/tests/test_loading.py | 8 ++-- sardem/utils.py | 9 ++++ setup.py | 2 +- 8 files changed, 101 insertions(+), 70 deletions(-) create mode 100644 sardem/constants.py diff --git a/sardem/cli.py b/sardem/cli.py index 254bb2e..b7ce0f7 100644 --- a/sardem/cli.py +++ b/sardem/cli.py @@ -1,15 +1,12 @@ """ Command line interface for running sardem """ -from sardem.download import Downloader import json -from argparse import ( - ArgumentError, - ArgumentParser, - ArgumentTypeError, - FileType, - RawTextHelpFormatter, -) +from argparse import (ArgumentError, ArgumentParser, ArgumentTypeError, + FileType, RawTextHelpFormatter) + +from sardem.download import Downloader +from sardem import utils def positive_small_int(argstring): @@ -135,6 +132,14 @@ def get_cli_args(): " Default is GDAL's top-left edge convention." ), ) + parser.add_argument( + "--use-exact-bbox", + action="store_true", + help=( + "If --bbox is a set of integers, don't pad the bbox by half a pixel. " + "Otherwise, assumes the DEM should be made from whole tiles." + ), + ) return parser.parse_args() @@ -143,11 +148,11 @@ def cli(): args = get_cli_args() import sardem.dem - if args.left_lon and args.geojson or args.left_lon and args.bbox: + if (args.left_lon and args.geojson) or (args.left_lon and args.bbox): raise ArgumentError( args.geojson, - "Can only use one of positional arguments (left_lon top_lat dlon dlat) " - ", --geojson, or --bbox", + "Can only use one type of bounding box specifier: (left_lon top_lat dlon dlat) " + ", --geojson, --bbox, or --wkt-file", ) # Need all 4 positionals, or the --geosjon elif ( @@ -156,19 +161,15 @@ def cli(): and not args.bbox and not args.wkt_file ): - raise ValueError("Need --bbox, --geojoin, or --wkt-file") + raise ValueError("Need --bbox, --geojson, or --wkt-file") geojson_dict = json.load(args.geojson) if args.geojson else None - if args.bbox: - left, bot, right, top = args.bbox - left_lon, top_lat = left, top - dlon = right - left - dlat = top - bot - elif args.wkt_file: - left_lon, top_lat, dlon, dlat = None, None, None, None - elif args.left_lon: + if args.left_lon: left_lon, top_lat = args.left_lon, args.top_lat dlon, dlat = args.dlon, args.dlat + bbox = utils.bounding_box(left_lon, top_lat, dlon, dlat) + elif args.bbox: + bbox = args.bbox if not args.output: output = ( @@ -178,10 +179,7 @@ def cli(): output = args.output sardem.dem.main( - left_lon, - top_lat, - dlon, - dlat, + bbox=bbox, geojson=geojson_dict, wkt_file=args.wkt_file, data_source=args.data_source, @@ -190,5 +188,6 @@ def cli(): make_isce_xml=args.make_isce_xml, keep_egm=args.keep_egm, shift_rsc=args.shift_rsc, + use_exact_bbox=args.use_exact_bbox, output_name=output, ) diff --git a/sardem/constants.py b/sardem/constants.py new file mode 100644 index 0000000..00e7a50 --- /dev/null +++ b/sardem/constants.py @@ -0,0 +1,2 @@ +NUM_PIXELS_SRTM1 = 3601 # For SRTM1 +DEFAULT_RES = 1 / 3600.0 diff --git a/sardem/cop_dem.py b/sardem/cop_dem.py index 40de034..dc5cee5 100644 --- a/sardem/cop_dem.py +++ b/sardem/cop_dem.py @@ -4,10 +4,11 @@ import requests from sardem import conversions, utils +from sardem.constants import DEFAULT_RES 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) @@ -44,6 +45,7 @@ def download_and_stitch( t_srs = 'epsg:4326' xres = DEFAULT_RES / xrate yres = DEFAULT_RES / yrate + resamp = "bilinear" if (xrate > 1 or yrate > 1) else "nearest" # access_mode = "overwrite" if overwrite else None option_dict = dict( @@ -54,7 +56,7 @@ def download_and_stitch( xRes=xres, yRes=yres, outputType=gdal.GDT_Int16, - resampleAlg="bilinear", + resampleAlg=resamp, multithread=True, warpMemoryLimit=5000, warpOptions=["NUM_THREADS=4"], @@ -82,7 +84,7 @@ def _gdal_cmd_from_options(src, dst, option_dict): from osgeo import gdal opts = deepcopy(option_dict) - # To see what the list of cli options are + # To see what the list of cli options are (gdal >= 3.5.0) opts["options"] = "__RETURN_OPTION_LIST__" opt_list = gdal.WarpOptions(**opts) out_opt_list = deepcopy(opt_list) diff --git a/sardem/dem.py b/sardem/dem.py index 0752770..c575257 100644 --- a/sardem/dem.py +++ b/sardem/dem.py @@ -1,4 +1,4 @@ -"""Digital Elevation Map (DEM) downloading/stitching/upsampling +"""Digital Elevation Map (DEM) downloading, stitching, upsampling Module contains utilities for downloading all necessary .hgt files for a lon/lat rectangle, stiches them into one DEM, and creates a @@ -38,15 +38,17 @@ PROJECTION LL """ from __future__ import division, print_function + import collections import logging import os + import numpy as np -from sardem import utils, loading, upsample_cy, conversions -from sardem.download import Tile, Downloader +from sardem import conversions, loading, upsample_cy, utils +from sardem.constants import NUM_PIXELS_SRTM1 +from sardem.download import Downloader, Tile -NUM_PIXELS = 3601 # For SRTM1 RSC_KEYS = [ "WIDTH", "FILE_LENGTH", @@ -79,7 +81,7 @@ class Stitcher: """ def __init__( - self, tile_names, filenames=[], data_source="NASA", num_pixels=NUM_PIXELS + self, tile_names, filenames=[], data_source="NASA", num_pixels=NUM_PIXELS_SRTM1 ): """List should come from Tile.srtm1_tile_names()""" self.tile_file_list = list(tile_names) @@ -191,7 +193,7 @@ def _load_tile(self, tile_name): else: return loading.load_elevation(filename) else: - return np.zeros((NUM_PIXELS, NUM_PIXELS), dtype=self.dtype) + return np.zeros((NUM_PIXELS_SRTM1, NUM_PIXELS_SRTM1), dtype=self.dtype) def load_and_stitch(self): """Function to load combine .hgt tiles @@ -279,11 +281,11 @@ def create_dem_rsc(self): return rsc_dict -def crop_stitched_dem(bounds, stitched_dem, rsc_data): - """Takes the output of Stitcher.load_and_stitch, crops to bounds +def crop_stitched_dem(bbox, stitched_dem, rsc_data): + """Takes the output of Stitcher.load_and_stitch, crops to bbox Args: - bounds (tuple[float]): (left, bot, right, top) lats and lons of + bbox (tuple[float]): (left, bot, right, top) lats and lons of desired bounding box for the DEM stitched_dem (numpy.array, 2D): result from files through Stitcher.load_and_stitch() @@ -293,7 +295,7 @@ def crop_stitched_dem(bounds, stitched_dem, rsc_data): numpy.array: a cropped version of the bigger stitched_dem """ indexes, new_starts = utils.find_bounding_idxs( - bounds, + bbox, rsc_data["X_STEP"], rsc_data["Y_STEP"], rsc_data["X_FIRST"], @@ -305,11 +307,12 @@ def crop_stitched_dem(bounds, stitched_dem, rsc_data): return cropped_dem, new_starts, new_sizes +def _float_is_on_bounds(x): + return int(x) == x + + def main( - left_lon=None, - top_lat=None, - dlon=None, - dlat=None, + bbox=None, geojson=None, wkt_file=None, data_source=None, @@ -318,16 +321,16 @@ def main( make_isce_xml=False, keep_egm=False, shift_rsc=False, + use_exact_bbox=False, output_name=None, ): """Function for entry point to create a DEM with `sardem` Args: - left_lon (float): Left most longitude of DEM box - top_lat (float): Top most longitude of DEM box - dlon (float): Width of box in longitude degrees - dlat (float): Height of box in latitude degrees - geojson (dict): geojson object outlining DEM (alternative to lat/lon) + bbox (tuple[float]): (left, bot, right, top) + Longitude/latitude desired bounding box for the DEM + geojson (dict): geojson object outlining DEM (alternative to bbox) + wkt_file (str): path to .wkt file outlining DEM (alternative to bbox) data_source (str): 'NASA' or 'AWS', where to download .hgt tiles from xrate (int): x-rate (columns) to upsample DEM (positive int) yrate (int): y-rate (rows) to upsample DEM (positive int) @@ -337,21 +340,33 @@ def main( shift_rsc (bool): Shift the .dem.rsc file down/right so that the X_FIRST and Y_FIRST values represent the pixel *center* (instead of GDAL's convention of pixel edge). Default = False. + use_exact_bbox (bool): If the bbox is a set of integers, avoid padding so that + whole tiles are used. + Default = False, so the bbox is padded to the nearest tile. output_name (str): name of file to save final DEM (default = elevation.dem) """ - if geojson: - bounds = utils.bounding_box(geojson=geojson) - elif wkt_file: - bounds = utils.get_wkt_bbox(wkt_file) - else: - bounds = utils.bounding_box(left_lon, top_lat, dlon, dlat) - logger.info("Bounds: %s", " ".join(str(b) for b in bounds)) - outrows, outcols = utils.get_output_size(bounds, xrate, yrate) + if bbox is None: + if geojson: + bbox = utils.bounding_box(geojson=geojson) + elif wkt_file: + bbox = utils.get_wkt_bbox(wkt_file) + + if bbox is None: + raise ValueError("Must provide either bbox or geojson or wkt_file") + logger.info("Bounds: %s", " ".join(str(b) for b in bbox)) + + if not use_exact_bbox: + if all(_float_is_on_bounds(b) for b in bbox): + logger.info("Shifting bbox to nearest tile bounds") + bbox = utils.shift_integer_bbox(bbox) + logger.info("New edge bounds: %s", " ".join(str(b) for b in bbox)) + + outrows, outcols = utils.get_output_size(bbox, xrate, yrate) if outrows * outcols > WARN_LIMIT: logger.warning( "Caution: Output size is {} x {} pixels.".format(outrows, outcols) ) - logger.warning("Are the bounds correct?") + logger.warning("Are the bounds correct (left, bottom, right, top)?") # Are we using GDAL's convention (pixel edge) or the center? # i.e. if `shift_rsc` is False, then we are `using_gdal_bounds` @@ -363,7 +378,7 @@ def main( cop_dem.download_and_stitch( output_name, - bounds, + bbox, keep_egm=keep_egm, xrate=xrate, yrate=yrate, @@ -375,7 +390,7 @@ def main( ) return - tile_names = list(Tile(*bounds).srtm1_tile_names()) + tile_names = list(Tile(*bbox).srtm1_tile_names()) d = Downloader(tile_names, data_source=data_source) local_filenames = d.download_all() @@ -386,10 +401,10 @@ def main( # Now create corresponding rsc file rsc_dict = s.create_dem_rsc() - # Cropping: get very close to the bounds asked for: + # Cropping: get very close to the bbox asked for: logger.info("Cropping stitched DEM to boundaries") stitched_dem, new_starts, new_sizes = crop_stitched_dem( - bounds, stitched_dem, rsc_dict + bbox, stitched_dem, rsc_dict ) new_x_first, new_y_first = new_starts new_rows, new_cols = new_sizes diff --git a/sardem/tests/test_dem.py b/sardem/tests/test_dem.py index aa8e8e7..d413aaf 100644 --- a/sardem/tests/test_dem.py +++ b/sardem/tests/test_dem.py @@ -1,12 +1,12 @@ -import unittest -import json -import tempfile -import shutil -from os.path import join, dirname import os +import shutil +import tempfile +import unittest +from os.path import dirname, join + import responses -from sardem import dem, utils, download +from sardem import dem, download, utils DATAPATH = join(dirname(__file__), "data") NETRC_PATH = join(DATAPATH, "netrc") @@ -172,3 +172,7 @@ def test_bounding_box(self): output_dem = sario.load_file('elevation.dem') # assert_array_almost_equal(expected_dem) """ + + +# TODO: tests with the hawaii COP tile +# https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N19_00_W156_00_DEM/Copernicus_DSM_COG_10_N19_00_W156_00_DEM.tif \ No newline at end of file diff --git a/sardem/tests/test_loading.py b/sardem/tests/test_loading.py index ca4dc75..dcd2375 100644 --- a/sardem/tests/test_loading.py +++ b/sardem/tests/test_loading.py @@ -1,9 +1,9 @@ -from collections import OrderedDict -import unittest -from os.path import join, dirname -import tempfile import shutil +import tempfile +import unittest import zipfile +from collections import OrderedDict +from os.path import dirname, join from sardem import loading diff --git a/sardem/utils.py b/sardem/utils.py index 6766684..a9e34f6 100644 --- a/sardem/utils.py +++ b/sardem/utils.py @@ -8,6 +8,7 @@ from math import floor from sardem import loading +from sardem.constants import DEFAULT_RES def set_logger_handler(logger, level="INFO"): @@ -112,6 +113,14 @@ def bounding_box(left_lon=None, top_lat=None, dlon=None, dlat=None, geojson=None return left, bottom, right, top +def shift_integer_bbox(bbox): + """Shift the integer bounds of a bbox by 1/2 pixel to select a whole tile""" + left, bottom, right, top = bbox + hp = 0.5 * DEFAULT_RES # half pixel + # Tile names refer to the center of the bottom-left corner of the tile + return left - hp, bottom + hp, right - hp, top + hp + + def coords(geojson): """Finds the coordinates of a geojson polygon Note: we are assuming one simple polygon with no holes diff --git a/setup.py b/setup.py index 3722ffe..3243751 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="sardem", - version="0.10.7", + version="0.10.8", author="Scott Staniewicz", author_email="scott.stanie@gmail.com", description="Create upsampled DEMs for InSAR processing", From 0877acf9dbcbbb3fe75aba08198f1e75c30bf0b7 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 10:50:32 -0400 Subject: [PATCH 02/11] round isce xml coordinates to 9 decimals --- sardem/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sardem/utils.py b/sardem/utils.py index a9e34f6..d3b858b 100644 --- a/sardem/utils.py +++ b/sardem/utils.py @@ -442,8 +442,8 @@ def gdal2isce_xml(fname, keep_egm=False, using_gdal_bounds=True): first_lon += 0.5 * delta_lon first_lat += 0.5 * delta_lat - img.firstLongitude = first_lon - img.firstLatitude = first_lat + img.firstLongitude = round(first_lon, 9) # rounding to avoid precision issues + img.firstLatitude = round(first_lat, 9) img.deltaLongitude = delta_lon img.deltaLatitude = delta_lat From c3ae4355c6f8b3f65a21aca89c26a3268bf4144f Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 11:11:17 -0400 Subject: [PATCH 03/11] use the original integer bbox to download otherwise, 3 extra ones are downloaded that arent used --- sardem/dem.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/sardem/dem.py b/sardem/dem.py index c575257..68ce6fd 100644 --- a/sardem/dem.py +++ b/sardem/dem.py @@ -355,6 +355,7 @@ def main( raise ValueError("Must provide either bbox or geojson or wkt_file") logger.info("Bounds: %s", " ".join(str(b) for b in bbox)) + bbox_orig = bbox if not use_exact_bbox: if all(_float_is_on_bounds(b) for b in bbox): logger.info("Shifting bbox to nearest tile bounds") @@ -372,6 +373,7 @@ def main( # i.e. if `shift_rsc` is False, then we are `using_gdal_bounds` using_gdal_bounds = not shift_rsc + # For copernicus, use GDAL to warp from the VRT if data_source == "COP": utils._gdal_installed_correctly() from sardem import cop_dem @@ -390,7 +392,8 @@ def main( ) return - tile_names = list(Tile(*bbox).srtm1_tile_names()) + # If using SRTM, download tiles manually and stitch + tile_names = list(Tile(*bbox_orig).srtm1_tile_names()) d = Downloader(tile_names, data_source=data_source) local_filenames = d.download_all() From 15b02ec2b2fd606dc58dd273e32884bb4cf6a493 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 21:18:40 -0400 Subject: [PATCH 04/11] use bbox orig for cropping --- sardem/dem.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/sardem/dem.py b/sardem/dem.py index 68ce6fd..f0d01bf 100644 --- a/sardem/dem.py +++ b/sardem/dem.py @@ -350,11 +350,11 @@ def main( bbox = utils.bounding_box(geojson=geojson) elif wkt_file: bbox = utils.get_wkt_bbox(wkt_file) - + if bbox is None: raise ValueError("Must provide either bbox or geojson or wkt_file") logger.info("Bounds: %s", " ".join(str(b) for b in bbox)) - + bbox_orig = bbox if not use_exact_bbox: if all(_float_is_on_bounds(b) for b in bbox): @@ -407,7 +407,7 @@ def main( # Cropping: get very close to the bbox asked for: logger.info("Cropping stitched DEM to boundaries") stitched_dem, new_starts, new_sizes = crop_stitched_dem( - bbox, stitched_dem, rsc_dict + bbox_orig, stitched_dem, rsc_dict ) new_x_first, new_y_first = new_starts new_rows, new_cols = new_sizes @@ -416,8 +416,13 @@ def main( rsc_dict["Y_FIRST"] = new_y_first rsc_dict["FILE_LENGTH"] = new_rows rsc_dict["WIDTH"] = new_cols - if shift_rsc: - rsc_dict = utils.shift_rsc_dict(rsc_dict, to_gdal=True) + + print( + f"{bbox = }, {bbox_orig = }, {rsc_dict['Y_FIRST'] = }, {rsc_dict['X_FIRST'] = }" + ) + + if not using_gdal_bounds: + rsc_dict = utils.shift_rsc_dict(rsc_dict, to_gdal=False) rsc_filename = output_name + ".rsc" From b784073459ff72099cbf999bf16c6de1ae4207de Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 22:21:40 -0400 Subject: [PATCH 05/11] fix shifting of bbox to happen onces at end --- sardem/conversions.py | 10 ++------ sardem/dem.py | 37 ++++++++++----------------- sardem/download.py | 12 ++++++--- sardem/utils.py | 59 ++++++++++++++++++++++++++++++++----------- 4 files changed, 69 insertions(+), 49 deletions(-) diff --git a/sardem/conversions.py b/sardem/conversions.py index 1c7920d..ae5ec08 100644 --- a/sardem/conversions.py +++ b/sardem/conversions.py @@ -68,7 +68,7 @@ def _get_size(filename): return xsize, ysize -def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_bounds=False): +def convert_dem_to_wgs84(dem_filename, geoid="egm96"): """Convert the file `dem_filename` from EGM96 heights to WGS84 ellipsoidal heights Overwrites file, requires GDAL to be installed @@ -79,8 +79,6 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_bounds=False): path_, fname = os.path.split(dem_filename) rsc_filename = os.path.join(path_, fname + ".rsc") - if not using_gdal_bounds: - utils.shift_rsc_file(rsc_filename, to_gdal=True) output_egm = os.path.join(path_, "egm_" + fname) # output_wgs = dem_filename.replace(ext, ".wgs84" + ext) @@ -97,8 +95,4 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_bounds=False): logger.error("Failed to convert DEM:", exc_info=True) logger.error("Reverting back, using EGM dem as output") os.rename(output_egm, dem_filename) - os.rename(rsc_filename_egm, rsc_filename) - - 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) + os.rename(rsc_filename_egm, rsc_filename) \ No newline at end of file diff --git a/sardem/dem.py b/sardem/dem.py index f0d01bf..6b77ebd 100644 --- a/sardem/dem.py +++ b/sardem/dem.py @@ -46,7 +46,7 @@ import numpy as np from sardem import conversions, loading, upsample_cy, utils -from sardem.constants import NUM_PIXELS_SRTM1 +from sardem.constants import DEFAULT_RES, NUM_PIXELS_SRTM1 from sardem.download import Downloader, Tile RSC_KEYS = [ @@ -355,13 +355,15 @@ def main( raise ValueError("Must provide either bbox or geojson or wkt_file") logger.info("Bounds: %s", " ".join(str(b) for b in bbox)) - bbox_orig = bbox if not use_exact_bbox: if all(_float_is_on_bounds(b) for b in bbox): logger.info("Shifting bbox to nearest tile bounds") bbox = utils.shift_integer_bbox(bbox) logger.info("New edge bounds: %s", " ".join(str(b) for b in bbox)) + + # Now we're assuming that `bbox` refers to the edges of the desired bounding box + # Print a warning if they're possibly requesting too-large a box by mistake outrows, outcols = utils.get_output_size(bbox, xrate, yrate) if outrows * outcols > WARN_LIMIT: logger.warning( @@ -369,10 +371,6 @@ def main( ) logger.warning("Are the bounds correct (left, bottom, right, top)?") - # Are we using GDAL's convention (pixel edge) or the center? - # i.e. if `shift_rsc` is False, then we are `using_gdal_bounds` - using_gdal_bounds = not shift_rsc - # For copernicus, use GDAL to warp from the VRT if data_source == "COP": utils._gdal_installed_correctly() @@ -387,13 +385,11 @@ def main( ) if make_isce_xml: logger.info("Creating ISCE2 XML file") - utils.gdal2isce_xml( - output_name, keep_egm=keep_egm, using_gdal_bounds=using_gdal_bounds - ) + utils.gdal2isce_xml(output_name, keep_egm=keep_egm) return # If using SRTM, download tiles manually and stitch - tile_names = list(Tile(*bbox_orig).srtm1_tile_names()) + tile_names = list(Tile(*bbox).srtm1_tile_names()) d = Downloader(tile_names, data_source=data_source) local_filenames = d.download_all() @@ -407,7 +403,7 @@ def main( # Cropping: get very close to the bbox asked for: logger.info("Cropping stitched DEM to boundaries") stitched_dem, new_starts, new_sizes = crop_stitched_dem( - bbox_orig, stitched_dem, rsc_dict + bbox, stitched_dem, rsc_dict ) new_x_first, new_y_first = new_starts new_rows, new_cols = new_sizes @@ -417,12 +413,7 @@ def main( rsc_dict["FILE_LENGTH"] = new_rows rsc_dict["WIDTH"] = new_cols - print( - f"{bbox = }, {bbox_orig = }, {rsc_dict['Y_FIRST'] = }, {rsc_dict['X_FIRST'] = }" - ) - - if not using_gdal_bounds: - rsc_dict = utils.shift_rsc_dict(rsc_dict, to_gdal=False) + rsc_dict = utils.shift_rsc_dict(rsc_dict, to_gdal=True) rsc_filename = output_name + ".rsc" @@ -475,17 +466,17 @@ def main( if make_isce_xml: logger.info("Creating ISCE2 XML file") - utils.gdal2isce_xml( - output_name, keep_egm=keep_egm, using_gdal_bounds=using_gdal_bounds - ) + utils.gdal2isce_xml(output_name, keep_egm=keep_egm) if keep_egm or data_source == "NASA_WATER": logger.info("Keeping DEM as EGM96 geoid heights") else: logger.info("Correcting DEM to heights above WGS84 ellipsoid") - conversions.convert_dem_to_wgs84( - output_name, using_gdal_bounds=using_gdal_bounds - ) + conversions.convert_dem_to_wgs84(output_name, geoid="egm96") + + # If the user wants the .rsc file to point to pixel center: + if shift_rsc: + utils.shift_rsc_file(rsc_filename, to_gdal=False) # Overwrite with smaller dtype for water mask if data_source == "NASA_WATER": diff --git a/sardem/download.py b/sardem/download.py index 7e506e3..422a45a 100644 --- a/sardem/download.py +++ b/sardem/download.py @@ -10,7 +10,8 @@ import numpy as np import requests -from . import utils +from sardem import utils +from sardem.constants import DEFAULT_RES try: input = raw_input # Check for python 2 @@ -145,11 +146,16 @@ def srtm1_tile_names(self): True >>> list(d.srtm1_tile_names()) ['N19W156', 'N19W155'] - >>> bounds = [-156.0, 19.0, -154.0, 20.0] # Show int bounds + >>> hp = 1 / 3600 / 2 # Half a degree in decimal degrees + >>> bounds = [-156.0 - hp, 19.0 - hp, -154.0 + hp, 20.0 + hp] >>> list(d.srtm1_tile_names()) ['N19W156', 'N19W155'] """ - left, bottom, right, top = self.bounds + # `bounds` should refer to the edges of the bounding box + # shift each inward by half pixel so they point to the boundary pixel centers + shift_inward = np.array([1, 1, -1, -1]) * DEFAULT_RES / 2 + left, bottom, right, top = np.array(self.bounds) + shift_inward + left_int, top_int = self.srtm1_tile_corner(left, top) right_int, bot_int = self.srtm1_tile_corner(right, bottom) # If exact integer was requested for top/right, assume tile with that number diff --git a/sardem/utils.py b/sardem/utils.py index d3b858b..c4e8585 100644 --- a/sardem/utils.py +++ b/sardem/utils.py @@ -145,36 +145,62 @@ def coords(geojson): return geojson["coordinates"][0] -def find_bounding_idxs(bounds, x_step, y_step, x_first, y_first): +def find_bounding_idxs(bbox, x_step, y_step, x_first, y_first): """Finds the indices of stitched dem to crop bounding box Also finds the new x_start and y_start after cropping. - Note: x_start, y_start could be different from bounds + Note: x_start, y_start could be different from bbox if steps didnt exactly match, but should be further up and left - Takes the desired bounds, .rsc data from stitched dem, + Args: + bbox (tuple[float]): (left,bottom,right,top) of bounding box + This refers to the *edges* of the box + x_step (float): step size in x direction + y_step (float): step size in y direction + x_first (float): x coordinate of first point in x direction + This refers to the *center* of the first pixel + y_first (float): y coordinate of first point in y direction + This refers to the *center* of the first pixel + + Takes the desired bbox, .rsc data from stitched dem, Examples: - >>> bounds = [-155.49, 19.0, -154.5, 19.5] + >>> bbox = [-155.05, 19.05, -154.05, 20.05] >>> x_step = 0.1 >>> y_step = -0.1 - >>> x_first = -156 + >>> x_first = -155 >>> y_first = 20.0 - >>> print(find_bounding_idxs(bounds, x_step, y_step, x_first, y_first)) - ((5, 10, 15, 5), (-155.5, 19.5)) - >>> bounds[-1] += 0.06 - >>> print(find_bounding_idxs(bounds, x_step, y_step, x_first, y_first)) - ((5, 10, 15, 4), (-155.5, 19.6)) + >>> print(find_bounding_idxs(bbox, x_step, y_step, x_first, y_first)) + ((0, 10, 10, 0), (-155.0, 20.0)) + >>> bbox[-1] -= 0.1 + >>> print(find_bounding_idxs(bbox, x_step, y_step, x_first, y_first)) + ((0, 10, 10, 1), (-155.0, 19.9)) """ + # `bbox` should refer to the edges of the bounding box + # shift by half pixel so they point to the pixel centers for index finding + hp = 0.5 * DEFAULT_RES # half pixel + left, bot, right, top = bbox + left += hp + bot += hp + + # Shift these two inward to be the final pixel centers + right -= hp + top -= hp - left, bot, right, top = bounds left_idx = int(round((left - x_first) / x_step)) - right_idx = int(round((right - x_first) / x_step)) + new_x_first = x_first + x_step * left_idx + right_idx = int(round((right - x_first) / x_step)) + 1 + # Note: y_step will be negative for these top_idx = int(round((top - y_first) / y_step)) - bot_idx = int(round((bot - y_first) / y_step)) - new_x_first = x_first + x_step * left_idx new_y_first = y_first + y_step * top_idx # Again: y_step negative + bot_idx = int(round((bot - y_first) / y_step)) + 1 + if any(arg < 0 for arg in (left_idx, right_idx, top_idx, bot_idx)): + raise ValueError( + "x_first/y_first ({}, {}) must be within the bbox {}".format( + x_first, y_first, bbox + ) + ) return (left_idx, bot_idx, right_idx, top_idx), (new_x_first, new_y_first) @@ -357,6 +383,7 @@ def gdal2isce_xml(fname, keep_egm=False, using_gdal_bounds=True): """ _gdal_installed_correctly() from osgeo import gdal + try: import isce # noqa import isceobj @@ -463,7 +490,9 @@ def _add_reference_datum(xml_file, keep_egm=False): import xml.etree.ElementTree as ET from xml.dom import minidom - logger.info("add info to xml file: {}".format(os.path.basename(xml_file))) + logger.info( + "add info to xml file: {}".format(os.path.basename(xml_file)) + ) # get property element for reference ref = ET.Element("property", attrib={"name": "reference"}) From 87014fe6552f2e945f73450aba80fefb2fdbf6f0 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 22:29:00 -0400 Subject: [PATCH 06/11] remove exact bbox cli arg not sure who would ever want that --- sardem/cli.py | 10 ---------- sardem/dem.py | 16 +++++----------- 2 files changed, 5 insertions(+), 21 deletions(-) diff --git a/sardem/cli.py b/sardem/cli.py index b7ce0f7..1810e8b 100644 --- a/sardem/cli.py +++ b/sardem/cli.py @@ -132,15 +132,6 @@ def get_cli_args(): " Default is GDAL's top-left edge convention." ), ) - parser.add_argument( - "--use-exact-bbox", - action="store_true", - help=( - "If --bbox is a set of integers, don't pad the bbox by half a pixel. " - "Otherwise, assumes the DEM should be made from whole tiles." - ), - ) - return parser.parse_args() @@ -188,6 +179,5 @@ def cli(): make_isce_xml=args.make_isce_xml, keep_egm=args.keep_egm, shift_rsc=args.shift_rsc, - use_exact_bbox=args.use_exact_bbox, output_name=output, ) diff --git a/sardem/dem.py b/sardem/dem.py index 6b77ebd..0a8b0a4 100644 --- a/sardem/dem.py +++ b/sardem/dem.py @@ -321,7 +321,6 @@ def main( make_isce_xml=False, keep_egm=False, shift_rsc=False, - use_exact_bbox=False, output_name=None, ): """Function for entry point to create a DEM with `sardem` @@ -340,9 +339,6 @@ def main( shift_rsc (bool): Shift the .dem.rsc file down/right so that the X_FIRST and Y_FIRST values represent the pixel *center* (instead of GDAL's convention of pixel edge). Default = False. - use_exact_bbox (bool): If the bbox is a set of integers, avoid padding so that - whole tiles are used. - Default = False, so the bbox is padded to the nearest tile. output_name (str): name of file to save final DEM (default = elevation.dem) """ if bbox is None: @@ -355,12 +351,10 @@ def main( raise ValueError("Must provide either bbox or geojson or wkt_file") logger.info("Bounds: %s", " ".join(str(b) for b in bbox)) - if not use_exact_bbox: - if all(_float_is_on_bounds(b) for b in bbox): - logger.info("Shifting bbox to nearest tile bounds") - bbox = utils.shift_integer_bbox(bbox) - logger.info("New edge bounds: %s", " ".join(str(b) for b in bbox)) - + if all(_float_is_on_bounds(b) for b in bbox): + logger.info("Shifting bbox to nearest tile bounds") + bbox = utils.shift_integer_bbox(bbox) + logger.info("New edge bounds: %s", " ".join(str(b) for b in bbox)) # Now we're assuming that `bbox` refers to the edges of the desired bounding box # Print a warning if they're possibly requesting too-large a box by mistake @@ -391,7 +385,7 @@ def main( # If using SRTM, download tiles manually and stitch tile_names = list(Tile(*bbox).srtm1_tile_names()) - d = Downloader(tile_names, data_source=data_source) + d = Downloader(tile_names, data_source=data_source, cache_dir=cache_dir) local_filenames = d.download_all() s = Stitcher(tile_names, filenames=local_filenames, data_source=data_source) From d24ca132e1a0f880a7525c295ccc910748d1c072 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 22:33:30 -0400 Subject: [PATCH 07/11] add an option of where cache dir should be --- sardem/cli.py | 18 +++++++++++++++--- sardem/dem.py | 6 ++++-- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/sardem/cli.py b/sardem/cli.py index 1810e8b..270596c 100644 --- a/sardem/cli.py +++ b/sardem/cli.py @@ -2,8 +2,13 @@ Command line interface for running sardem """ import json -from argparse import (ArgumentError, ArgumentParser, ArgumentTypeError, - FileType, RawTextHelpFormatter) +from argparse import ( + ArgumentError, + ArgumentParser, + ArgumentTypeError, + FileType, + RawTextHelpFormatter, +) from sardem.download import Downloader from sardem import utils @@ -132,6 +137,12 @@ def get_cli_args(): " Default is GDAL's top-left edge convention." ), ) + parser.add_argument( + "--cache-dir", + help=( + "Location to save downloaded files (Default = {})".format(utils.get_cache_dir()) + ), + ) return parser.parse_args() @@ -170,6 +181,7 @@ def cli(): output = args.output sardem.dem.main( + output_name=output, bbox=bbox, geojson=geojson_dict, wkt_file=args.wkt_file, @@ -179,5 +191,5 @@ def cli(): make_isce_xml=args.make_isce_xml, keep_egm=args.keep_egm, shift_rsc=args.shift_rsc, - output_name=output, + cache_dir=args.cache_dir, ) diff --git a/sardem/dem.py b/sardem/dem.py index 0a8b0a4..3b574c6 100644 --- a/sardem/dem.py +++ b/sardem/dem.py @@ -312,6 +312,7 @@ def _float_is_on_bounds(x): def main( + output_name=None, bbox=None, geojson=None, wkt_file=None, @@ -321,11 +322,12 @@ def main( make_isce_xml=False, keep_egm=False, shift_rsc=False, - output_name=None, + cache_dir=None, ): """Function for entry point to create a DEM with `sardem` Args: + output_name (str): name of file to save final DEM (default = elevation.dem) bbox (tuple[float]): (left, bot, right, top) Longitude/latitude desired bounding box for the DEM geojson (dict): geojson object outlining DEM (alternative to bbox) @@ -339,7 +341,7 @@ def main( shift_rsc (bool): Shift the .dem.rsc file down/right so that the X_FIRST and Y_FIRST values represent the pixel *center* (instead of GDAL's convention of pixel edge). Default = False. - output_name (str): name of file to save final DEM (default = elevation.dem) + cache_dir (str): directory to cache downloaded tiles """ if bbox is None: if geojson: From a1231e184712c8b0cb9b381f836a45f978ffad3f Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 22:37:25 -0400 Subject: [PATCH 08/11] change some to pytest --- sardem/tests/test_dem.py | 59 ++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 33 deletions(-) diff --git a/sardem/tests/test_dem.py b/sardem/tests/test_dem.py index d413aaf..6d76fd4 100644 --- a/sardem/tests/test_dem.py +++ b/sardem/tests/test_dem.py @@ -3,6 +3,7 @@ import tempfile import unittest from os.path import dirname, join +import pytest import responses @@ -65,23 +66,18 @@ def test_download(self): self.assertTrue(os.path.exists(d._filepath(self.test_tile))) -class TestRsc(unittest.TestCase): - def setUp(self): - self.rsc_path = join(DATAPATH, "elevation.dem.rsc") +class TestRsc: + rsc_path = join(DATAPATH, "elevation.dem.rsc") def test_upsample_dem_rsc(self): # Test input checking - self.assertRaises( - ValueError, - utils.upsample_dem_rsc, - xrate=2, + with pytest.raises(ValueError): + utils.upsample_dem_rsc( xrate=2, rsc_dict={"something": 1}, rsc_filename=self.rsc_path, ) - self.assertRaises(ValueError, utils.upsample_dem_rsc, xrate=2) - self.assertRaises( - ValueError, utils.upsample_dem_rsc, rsc_filename=self.rsc_path - ) # Need rate + with pytest.raises(ValueError): + utils.upsample_dem_rsc(xrate=2) up_rsc = utils.upsample_dem_rsc(xrate=1, yrate=1, rsc_filename=self.rsc_path) expected = """\ @@ -126,36 +122,33 @@ def test_upsample_dem_rsc(self): Z_SCALE 1 PROJECTION LL """ - self.assertEqual(expected, up_rsc) - - -class TestBounds(unittest.TestCase): - def setUp(self): - self.coords = [ - [-156.0, 18.7], - [-154.6, 18.7], - [-154.6, 20.3], - [-156.0, 20.3], - [-156.0, 18.7], - ] - self.left_lon = -156.0 - self.top_lat = 20.3 - self.dlat = 1.6 - self.dlon = 1.4 + assert expected == up_rsc + + +class TestBounds: + coords = [ + [-156.0, 18.7], + [-154.6, 18.7], + [-154.6, 20.3], + [-156.0, 20.3], + [-156.0, 18.7], + ] + left_lon = -156.0 + top_lat = 20.3 + dlat = 1.6 + dlon = 1.4 def test_corner_input(self): result = utils.corner_coords(self.left_lon, self.top_lat, self.dlon, self.dlat) - self.assertEqual( - set(tuple(c) for c in result), set(tuple(c) for c in self.coords) - ) + assert set(tuple(c) for c in result) == set(tuple(c) for c in self.coords) def test_bounding_box(self): - self.assertEqual( - utils.bounding_box(self.left_lon, self.top_lat, self.dlon, self.dlat), - (-156.0, 18.7, -154.6, 20.3), + assert utils.bounding_box(self.left_lon, self.top_lat, self.dlon, self.dlat) == ( + (-156.0, 18.7, -154.6, 20.3) ) + """ TODO: finish cropped, upsampled dem, show this: From 908cb63f75f7e53e6980bbbcaa38d06cc1f109aa Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 22:58:21 -0400 Subject: [PATCH 09/11] start adding real dem test --- sardem/tests/test_dem.py | 55 +++++++++++++++++++++++++++++++++------- 1 file changed, 46 insertions(+), 9 deletions(-) diff --git a/sardem/tests/test_dem.py b/sardem/tests/test_dem.py index 6d76fd4..79616d8 100644 --- a/sardem/tests/test_dem.py +++ b/sardem/tests/test_dem.py @@ -1,4 +1,7 @@ import os +import numpy as np +import zipfile +import gzip import shutil import tempfile import unittest @@ -8,6 +11,8 @@ import responses from sardem import dem, download, utils +from sardem.constants import DEFAULT_RES +HALF_PIXEL = 0.5 * DEFAULT_RES DATAPATH = join(dirname(__file__), "data") NETRC_PATH = join(DATAPATH, "netrc") @@ -72,10 +77,11 @@ class TestRsc: def test_upsample_dem_rsc(self): # Test input checking with pytest.raises(ValueError): - utils.upsample_dem_rsc( xrate=2, - rsc_dict={"something": 1}, - rsc_filename=self.rsc_path, - ) + utils.upsample_dem_rsc( + xrate=2, + rsc_dict={"something": 1}, + rsc_filename=self.rsc_path, + ) with pytest.raises(ValueError): utils.upsample_dem_rsc(xrate=2) @@ -143,10 +149,9 @@ def test_corner_input(self): assert set(tuple(c) for c in result) == set(tuple(c) for c in self.coords) def test_bounding_box(self): - assert utils.bounding_box(self.left_lon, self.top_lat, self.dlon, self.dlat) == ( - (-156.0, 18.7, -154.6, 20.3) - ) - + assert utils.bounding_box( + self.left_lon, self.top_lat, self.dlon, self.dlat + ) == ((-156.0, 18.7, -154.6, 20.3)) """ @@ -167,5 +172,37 @@ def test_bounding_box(self): """ +@pytest.fixture +def sample_hgt_file(tmp_path): + tmpfile = tmp_path / "sample.hgt" + path = join(DATAPATH, "N19W156.hgt.zip") + with zipfile.ZipFile(path, "r") as zip_ref: + with open(tmpfile, "wb") as f: + f.write(zip_ref.read("N19W156.hgt")) + return tmpfile + + +def sample_cop_tile(tmp_path): + path = join(DATAPATH, "cop_tile_hawaii.dem.gz") + tmpfile = tmp_path / "cop_tile_hawaii.dem" + with gzip.open(path, "rb") as f_in: + with open(tmpfile, "wb") as f_out: + f_out.write(f_in.read()) + return tmpfile + # return np.fromfile(tmpfile, dtype=np.int16).reshape(3600, 3600) + + +class TestMain: + bbox = [-156.0, 19.0, -155.0, 20.0] + + def test_main(self, tmp_path, sample_hgt_file): + cache_dir = sample_hgt_file.parent + sample_hgt = np.fromfile(sample_hgt_file, dtype=np.int16).reshape(3601, 3601) + dem.main(output_name=tmp_path, bbox=self.bbox, keep_egm=True, data_source="NASA", cache_dir=cache_dir) + output = np.fromfile(tmp_path, dtype=np.int16).reshape(3600, 3600) + np.testing.assert_allclose(sample_hgt, output[:-1, :-1]) + + + # TODO: tests with the hawaii COP tile -# https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N19_00_W156_00_DEM/Copernicus_DSM_COG_10_N19_00_W156_00_DEM.tif \ No newline at end of file +# https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N19_00_W156_00_DEM/Copernicus_DSM_COG_10_N19_00_W156_00_DEM.tif From 882c229c0247cca6f5a18a13b2111e0f06270cd4 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 23:17:29 -0400 Subject: [PATCH 10/11] first srtm test working --- sardem/tests/test_dem.py | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/sardem/tests/test_dem.py b/sardem/tests/test_dem.py index 79616d8..34cdbd8 100644 --- a/sardem/tests/test_dem.py +++ b/sardem/tests/test_dem.py @@ -10,7 +10,7 @@ import responses -from sardem import dem, download, utils +from sardem import dem, download, utils, loading from sardem.constants import DEFAULT_RES HALF_PIXEL = 0.5 * DEFAULT_RES @@ -171,17 +171,6 @@ def test_bounding_box(self): # assert_array_almost_equal(expected_dem) """ - -@pytest.fixture -def sample_hgt_file(tmp_path): - tmpfile = tmp_path / "sample.hgt" - path = join(DATAPATH, "N19W156.hgt.zip") - with zipfile.ZipFile(path, "r") as zip_ref: - with open(tmpfile, "wb") as f: - f.write(zip_ref.read("N19W156.hgt")) - return tmpfile - - def sample_cop_tile(tmp_path): path = join(DATAPATH, "cop_tile_hawaii.dem.gz") tmpfile = tmp_path / "cop_tile_hawaii.dem" @@ -195,12 +184,20 @@ def sample_cop_tile(tmp_path): class TestMain: bbox = [-156.0, 19.0, -155.0, 20.0] - def test_main(self, tmp_path, sample_hgt_file): - cache_dir = sample_hgt_file.parent - sample_hgt = np.fromfile(sample_hgt_file, dtype=np.int16).reshape(3601, 3601) - dem.main(output_name=tmp_path, bbox=self.bbox, keep_egm=True, data_source="NASA", cache_dir=cache_dir) - output = np.fromfile(tmp_path, dtype=np.int16).reshape(3600, 3600) - np.testing.assert_allclose(sample_hgt, output[:-1, :-1]) + def test_main(self, tmp_path): + path = join(DATAPATH, "N19W156.hgt.zip") + # tmpfile = tmp_path / "N19W156.hgt.zip" + unzipfile = tmp_path / "N19W156.hgt" + with zipfile.ZipFile(path, "r") as zip_ref: + with open(unzipfile, "wb") as f: + f.write(zip_ref.read("N19W156.hgt")) + sample_hgt = loading.load_elevation(unzipfile) + sample_hgt[sample_hgt < -1000] = 0 + + tmp_output = tmp_path / "output.dem" + dem.main(output_name=str(tmp_output), bbox=self.bbox, keep_egm=True, data_source="NASA", cache_dir=str(tmp_path)) + output = np.fromfile(tmp_output, dtype=np.int16).reshape(3600, 3600) + np.testing.assert_allclose(sample_hgt[:-1, :-1], output) From 6d6e44dc9e514cefcd0210a78df1461f3a5a41b7 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sun, 14 Aug 2022 23:21:08 -0400 Subject: [PATCH 11/11] add cop test too format --- sardem/tests/test_dem.py | 64 ++++++++++++++++++---------------------- 1 file changed, 29 insertions(+), 35 deletions(-) diff --git a/sardem/tests/test_dem.py b/sardem/tests/test_dem.py index 34cdbd8..98309fe 100644 --- a/sardem/tests/test_dem.py +++ b/sardem/tests/test_dem.py @@ -12,6 +12,7 @@ from sardem import dem, download, utils, loading from sardem.constants import DEFAULT_RES + HALF_PIXEL = 0.5 * DEFAULT_RES DATAPATH = join(dirname(__file__), "data") @@ -154,52 +155,45 @@ def test_bounding_box(self): ) == ((-156.0, 18.7, -154.6, 20.3)) -""" -TODO: - finish cropped, upsampled dem, show this: - expected_dem = np.array( - [[2071, 2072, 2074, 2076, 2078, 2078, 2079, 2080, 2082], [ - 2071, 2072, 2073, 2075, 2076, 2077, 2078, 2079, 2081 - ], [2071, 2072, 2073, 2074, 2075, 2076, 2078, 2079, 2080], [ - 2071, 2071, 2072, 2073, 2074, 2075, 2077, 2078, 2080 - ], [2071, 2071, 2072, 2073, 2074, 2075, 2076, 2078, 2080], [ - 2071, 2071, 2072, 2072, 2073, 2074, 2076, 2077, 2079 - ], [2071, 2071, 2072, 2072, 2073, 2074, 2076, 2077, 2078]], - dtype='