adjust the isce XML file to use pixel centers

otherwise, the XML has the edges as "startingvalue", "xmax", etc

sardem/conversions.py

@@ -69,7 +69,7 @@ def _get_size(filename):
     return xsize, ysize
 
 
-def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_rsc=False):
+def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_bounds=False):
     """Convert the file `dem_filename` from EGM96 heights to WGS84 ellipsoidal heights
 
     Overwrites file, requires GDAL to be installed
@@ -80,7 +80,7 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_rsc=False):
 
     path_, fname = os.path.split(dem_filename)
     rsc_filename = os.path.join(path_, fname + ".rsc")
-    if not using_gdal_rsc:
+    if not using_gdal_bounds:
         utils.shift_rsc_file(rsc_filename, to_gdal=True)
 
     output_egm = os.path.join(path_, "egm_" + fname)
@@ -100,7 +100,7 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96", using_gdal_rsc=False):
         os.rename(output_egm, dem_filename)
         os.rename(rsc_filename_egm, rsc_filename)
 
-    if not using_gdal_rsc:
+    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)
 

sardem/dem.py

@@ -348,19 +348,31 @@ def main(
     logger.info("Bounds: %s", " ".join(str(b) for b in bounds))
     outrows, outcols = utils.get_output_size(bounds, xrate, yrate)
     if outrows * outcols > WARN_LIMIT:
-        logger.warning("Caution: Output size is {} x {} pixels.".format(outrows, outcols))
+        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
+
     if data_source == "COP":
         utils._gdal_installed_correctly()
         from sardem import cop_dem
 
         cop_dem.download_and_stitch(
-            output_name, bounds, keep_egm=keep_egm, xrate=xrate, yrate=yrate,
+            output_name,
+            bounds,
+            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)
+            utils.gdal2isce_xml(
+                output_name, keep_egm=keep_egm, using_gdal_bounds=using_gdal_bounds
+            )
         return
 
     tile_names = list(Tile(*bounds).srtm1_tile_names())
@@ -440,14 +452,15 @@ def main(
 
     if make_isce_xml:
         logger.info("Creating ISCE2 XML file")
-        utils.gdal2isce_xml(output_name, keep_egm=keep_egm)
+        utils.gdal2isce_xml(output_name, keep_egm=keep_egm, shift_rsc=using_gdal_bounds)
 
     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")
-        using_gdal_rsc = not shift_rsc
-        conversions.convert_dem_to_wgs84(output_name, using_gdal_rsc=using_gdal_rsc)
+        conversions.convert_dem_to_wgs84(
+            output_name, using_gdal_bounds=using_gdal_bounds
+        )
 
     # Overwrite with smaller dtype for water mask
     if data_source == "NASA_WATER":

sardem/utils.py

@@ -338,7 +338,7 @@ def _gdal_installed_correctly():
         return False
 
 
-def gdal2isce_xml(fname, keep_egm=False):
+def gdal2isce_xml(fname, keep_egm=False, using_gdal_bounds=True):
     """
     Generate ISCE xml file from gdal supported file
 
@@ -424,10 +424,19 @@ def gdal2isce_xml(fname, keep_egm=False):
             logger.info("Assuming default, BSQ")
             img.scheme = "BSQ"
 
-    img.firstLongitude = transform[0]
-    img.firstLatitude = transform[3]
-    img.deltaLatitude = transform[5]
-    img.deltaLongitude = transform[1]
+    first_lon = transform[0]
+    first_lat = transform[3]
+    delta_lat = transform[5]
+    delta_lon = transform[1]
+    if using_gdal_bounds:
+        # We need to shift the `first________` values to the middle of the top left pixel
+        first_lon += 0.5 * delta_lon
+        first_lat += 0.5 * delta_lat
+
+    img.firstLongitude = first_lon
+    img.firstLatitude = first_lat
+    img.deltaLongitude = delta_lon
+    img.deltaLatitude = delta_lat
 
     xml_file = outname + ".xml"
     img.dump(xml_file)

setup.py

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