Merge pull request #10 from scottstanie/bbox-integer

enhancement: make integer bbox align with tiles

Author: scottstanie

sardem/cli.py

@@ -1,7 +1,6 @@
 """
 Command line interface for running sardem
 """
-from sardem.download import Downloader
 import json
 from argparse import (
     ArgumentError,
@@ -11,6 +10,9 @@
     RawTextHelpFormatter,
 )
 
+from sardem.download import Downloader
+from sardem import utils
+
 
 def positive_small_int(argstring):
     try:
@@ -135,19 +137,24 @@ 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()
 
 
 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 +163,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 +181,8 @@ def cli():
         output = args.output
 
     sardem.dem.main(
-        left_lon,
-        top_lat,
-        dlon,
-        dlat,
+        output_name=output,
+        bbox=bbox,
         geojson=geojson_dict,
         wkt_file=args.wkt_file,
         data_source=args.data_source,
@@ -190,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,
     )

sardem/constants.py

@@ -0,0 +1,2 @@
+NUM_PIXELS_SRTM1 = 3601  # For SRTM1
+DEFAULT_RES = 1 / 3600.0

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)

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)

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 DEFAULT_RES, 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,13 @@ 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,
+    output_name=None,
+    bbox=None,
     geojson=None,
     wkt_file=None,
     data_source=None,
@@ -318,16 +322,16 @@ 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:
-        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)
+        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)
+        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,47 +341,53 @@ 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 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 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(
             "Caution: Output size is {} x {} pixels.".format(outrows, outcols)
         )
-        logger.warning("Are the bounds correct?")
-
-    # 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
+        logger.warning("Are the bounds correct (left, bottom, right, top)?")
 
+    # For copernicus, use GDAL to warp from the VRT
     if data_source == "COP":
         utils._gdal_installed_correctly()
         from sardem import cop_dem
 
         cop_dem.download_and_stitch(
             output_name,
-            bounds,
+            bbox,
             keep_egm=keep_egm,
             xrate=xrate,
             yrate=yrate,
         )
         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
 
-    tile_names = list(Tile(*bounds).srtm1_tile_names())
+    # 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)
@@ -386,10 +396,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
@@ -398,8 +408,8 @@ 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)
+
+    rsc_dict = utils.shift_rsc_dict(rsc_dict, to_gdal=True)
 
     rsc_filename = output_name + ".rsc"
 
@@ -452,17 +462,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":