Merge pull request #4 from scottstanie/geoid-default

Geoid default

data/elevation.hdr

@@ -0,0 +1,15 @@
+ENVI
+description = {
+elevation.dem}
+samples = 5401
+lines   = 5401
+bands   = 1
+header offset = 0
+file type = ENVI Standard
+data type = 2
+interleave = bsq
+byte order = 0
+map info = {Geographic Lat/Lon, 1, 1, -121.20000000224, 37.30000000196, 0.000277777776999999, 0.000277777776999999,WGS-84}
+coordinate system string = {GEOGCS["GCS_unknown",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]}
+band names = {
+Band 1}

data/hawaii-test/elevation_wgs84.hdr

@@ -0,0 +1,15 @@
+ENVI
+description = {
+elevation_wgs84.dem}
+samples = 5761
+lines   = 5761
+bands   = 1
+header offset = 0
+file type = ENVI Standard
+data type = 2
+interleave = bsq
+byte order = 0
+map info = {Geographic Lat/Lon, 1, 1, -156.300138890849, 20.3001388908485, 0.000277777777, 0.000277777777,WGS-84}
+coordinate system string = {GEOGCS["GCS_unknown",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]}
+band names = {
+Band 1}

sardem/__init__.py

@@ -1,3 +1,3 @@
-from . import dem
+# from . import dem
 from . import utils
 from . import loading

sardem/cli.py

@@ -1,7 +1,7 @@
 """
 Command line interface for running sardem
 """
-from sardem.dem import Downloader
+from sardem.download import Downloader
 import sys
 import json
 from argparse import (
@@ -11,7 +11,6 @@
     FileType,
     RawTextHelpFormatter,
 )
-import sardem
 
 
 def positive_small_int(argstring):
@@ -38,7 +37,7 @@ def positive_small_int(argstring):
     Also creates elevation.dem.rsc with start lat/lon, stride, and other info."""
 
 
-def cli():
+def get_cli_args():
     parser = ArgumentParser(
         prog="sardem",
         description=DESCRIPTION,
@@ -112,16 +111,21 @@ def cli():
         help="Source of SRTM data (default %(default)s). See README for more.",
     )
     parser.add_argument(
-        "--convert-to-wgs84",
-        "-c",
+        "--keep-egm96",
         action="store_true",
         help=(
-            "Convert the DEM heights from geoid heights above EGM96 "
-            "to heights above WGS84 ellipsoid"
+            "Keep the DEM heights as geoid heights above EGM96. "
+            "Default is to convert to WGS84 for InSAR processing."
         ),
     )
 
-    args = parser.parse_args()
+    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:
         raise ArgumentError(
             args.geojson,
@@ -135,8 +139,7 @@ def cli():
         and not args.bbox
         and not args.wkt_file
     ):
-        parser.print_usage(sys.stderr)
-        sys.exit(1)
+        raise ValueError("Need --bbox, --geojoin, or --wkt-file")
 
     geojson_dict = json.load(args.geojson) if args.geojson else None
     if args.bbox:
@@ -167,6 +170,6 @@ def cli():
         args.data_source,
         args.xrate,
         args.yrate,
-        args.convert_to_wgs84,
+        args.keep_egm96,
         output,
     )

sardem/conversions.py

@@ -7,18 +7,28 @@
 
 logger = logging.getLogger("sardem")
 
-EGM_FILE = os.path.join(utils.get_cache_dir(), "egm96_15.gtx")
-
-
-def egm96_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True):
+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,
+}
+EGM_FILES = {
+    "egm96": os.path.join(utils.get_cache_dir(), "egm96_15.gtx"),
+    "egm08": os.path.join(utils.get_cache_dir(), "egm08_25.gtx"),
+}
+
+
+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"""
 
     if output is None:
         ext = os.path.splitext(filename)[1]
         output = filename.replace(ext, ".wgs84" + ext)
 
-    if not os.path.exists(EGM_FILE):
-        download_egm96_grid()
+    egm_file = EGM_FILES[geoid]
+    if not os.path.exists(egm_file):
+        download_egm_grid(geoid=geoid)
 
     # Note: https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-tr
     # If not specified, gdalwarp will generate an output raster with xres=yres
@@ -34,7 +44,7 @@ def egm96_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True):
         overwrite="-overwrite" if overwrite else "",
         xres=xres,
         yres=yres,
-        egm_file=EGM_FILE,
+        egm_file=egm_file,
     )
     logger.info(cmd)
     subprocess.run(cmd, check=True, shell=True)
@@ -56,7 +66,7 @@ def _get_resolution(filename):
     return xres, yres
 
 
-def convert_dem_to_wgs84(dem_filename):
+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
@@ -75,7 +85,7 @@ def convert_dem_to_wgs84(dem_filename):
     os.rename(dem_filename, output_egm)
     os.rename(rsc_filename, rsc_filename_egm)
     try:
-        egm96_to_wgs84(output_egm, output=dem_filename, overwrite=True, copy_rsc=True)
+        egm_to_wgs84(output_egm, output=dem_filename, overwrite=True, copy_rsc=True, geoid=geoid)
         os.remove(output_egm)
         os.remove(rsc_filename_egm)
     except Exception:
@@ -88,14 +98,21 @@ def convert_dem_to_wgs84(dem_filename):
     shift_dem_rsc(rsc_filename, to_gdal=False)
 
 
-def download_egm96_grid():
-    url = "http://download.osgeo.org/proj/vdatum/egm96_15/egm96_15.gtx"
-    if os.path.exists(EGM_FILE):
-        logger.info("{} already exists, skipping.".format(EGM_FILE))
+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
 
-    logger.info("Downloading from {} to {}".format(url, EGM_FILE))
-    with open(EGM_FILE, "wb") as f:
+    logger.info("Downloading from {} to {}".format(url, egm_file))
+    with open(egm_file, "wb") as f:
         resp = requests.get(url)
         f.write(resp.content)
 

sardem/dem.py

@@ -324,7 +324,7 @@ def main(
     data_source=None,
     xrate=1,
     yrate=1,
-    convert_to_wgs84=False,
+    keep_egm96=False,
     output_name=None,
 ):
     """Function for entry point to create a DEM with `sardem`
@@ -338,8 +338,8 @@ def main(
         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)
-        convert_to_wgs84 (bool): Convert the DEM heights from geoid heights
-            above EGM96 to heights above WGS84 ellipsoid
+        keep_egm96 (bool): Don't convert the DEM heights from geoid heights
+            above EGM96 to heights above WGS84 ellipsoid (default = False)
         output_name (str): name of file to save final DEM (usually elevation.dem)
     """
     if geojson:
@@ -423,11 +423,11 @@ def main(
         os.remove(dem_filename_small)
         os.remove(rsc_filename_small)
 
-    if convert_to_wgs84:
+    if keep_egm96:
+        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)
-    else:
-        logger.info("Keeping DEM as EGM96 geoid heights")
 
     # Overwrite with smaller dtype for water mask
     if data_source == "NASA_WATER":

sardem/download.py

@@ -1,12 +1,15 @@
-import logging
 import getpass
+import logging
 import math
-from multiprocessing.pool import ThreadPool
 import netrc
 import os
 import re
 import subprocess
+from multiprocessing.pool import ThreadPool
+
+import numpy as np
 import requests
+
 from . import utils
 
 try:
@@ -351,8 +354,9 @@ def download_and_save(self, tile_name):
                     logger.warning("Cannot find url %s, using zeros for tile." % url)
                     # Raise only if we want to kill everything
                     # response.raise_for_status()
-                    # Return False so caller can track failed downloads
-                    return None
+                    local_filename = os.path.splitext(local_filename)[0]
+                    self._write_zeros(local_filename)
+                    return local_filename
 
                 f.write(response.content)
                 logger.info("Writing to {}".format(local_filename))
@@ -370,6 +374,9 @@ def _all_files_exist(self):
         filepaths = [self._filepath(tile_name) for tile_name in self.tile_names]
         return all(os.path.exists(f) for f in filepaths)
 
+    def _write_zeros(self, local_filename):
+        np.zeros((3601, 3601), dtype=np.int16).tofile(local_filename)
+
     def download_all(self):
         """Downloads and saves all tiles from tile list"""
         # Only need to get credentials for this case:

sardem/loading.py

@@ -1,10 +1,6 @@
 from __future__ import division
 import collections
 import os
-import numpy as np
-
-INT_16_LE = np.dtype("<i2")
-INT_16_BE = np.dtype(">i2")
 
 RSC_KEY_TYPES = [
     ("width", int),
@@ -45,6 +41,10 @@ def load_elevation(filename):
     so either manually set data(data == np.min(data)) = 0,
         data = np.clip(data, 0, None), or when plotting, plt.imshow(data, vmin=0)
     """
+    import numpy as np
+
+    INT_16_LE = np.dtype("<i2")
+    INT_16_BE = np.dtype(">i2")
     ext = os.path.splitext(filename)[1]
     data_type = INT_16_LE if ext == ".dem" else INT_16_BE
     data = np.fromfile(filename, dtype=data_type)
@@ -86,6 +86,7 @@ def load_watermask(filename):
     Reference:
     https://lpdaac.usgs.gov/products/srtmswbdv003/
     """
+    import numpy as np
     return np.fromfile(filename, dtype=np.uint8).reshape((3601, 3601))
 
 

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="sardem",
-    version="0.7.0",
+    version="0.8.0",
     author="Scott Staniewicz",
     author_email="scott.stanie@utexas.com",
     description="Create upsampled DEMs for InSAR processing",