Merge fix-dateline-clean: handle antimeridian-crossing bounding boxes (fix #18)

Author: scottstanie

sardem/cop_dem.py

@@ -32,19 +32,116 @@ def download_and_stitch(
         https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198
         https://copernicus-dem-30m.s3.amazonaws.com/readme.html
     """
+    import os
+    import tempfile
+
     from osgeo import gdal
 
     gdal.UseExceptions()
-    # TODO: does downloading make it run any faster?
-    # if download_vrt:
-    #     cache_dir = utils.get_cache_dir()
-    #     vrt_filename = os.path.join(cache_dir, "cop_global.vrt")
-    #     if not os.path.exists(vrt_filename):
-    #         make_cop_vrt(vrt_filename)
-    # else:
+
     if vrt_filename is None:
         vrt_filename = "/vsicurl/https://raw.githubusercontent.com/scottstanie/sardem/master/sardem/data/cop_global.vrt"  # noqa
 
+    bboxes = utils.check_dateline(bbox)
+
+    if len(bboxes) == 1:
+        _download_single_bbox(
+            output_name,
+            bbox,
+            vrt_filename,
+            keep_egm,
+            xrate,
+            yrate,
+            output_format,
+            output_type,
+        )
+        return
+
+    # Dateline crossing: download each sub-bbox, shift tiles, merge
+    logger.info(
+        "Dateline crossing detected, downloading {} separate regions".format(
+            len(bboxes)
+        )
+    )
+
+    temp_files = []
+    with tempfile.TemporaryDirectory() as tmpdir:
+        for idx, sub_bbox in enumerate(bboxes):
+            temp_file = os.path.join(tmpdir, "dem_part_{}.tif".format(idx))
+            temp_files.append(temp_file)
+            logger.info("Downloading region {} of {}".format(idx + 1, len(bboxes)))
+            _download_single_bbox(
+                temp_file,
+                sub_bbox,
+                vrt_filename,
+                keep_egm,
+                xrate,
+                yrate,
+                "GTiff",
+                output_type,
+            )
+
+        # Shift eastern tiles so they're adjacent to western tiles in pixel space
+        for temp_file in temp_files:
+            _shift_tile_if_needed(temp_file)
+
+        logger.info("Merging {} regions into final DEM".format(len(temp_files)))
+        vrt_temp = os.path.join(tmpdir, "merged.vrt")
+        gdal.BuildVRT(vrt_temp, temp_files)
+
+        if output_format == "GTiff":
+            gdal.Warp(
+                output_name,
+                vrt_temp,
+                options=gdal.WarpOptions(
+                    format=output_format,
+                    multithread=True,
+                    callback=gdal.TermProgress,
+                ),
+            )
+        else:
+            gdal.Translate(
+                output_name,
+                vrt_temp,
+                format=output_format,
+                callback=gdal.TermProgress,
+            )
+
+
+def _shift_tile_if_needed(filepath):
+    """Shift tile geotransform by -360 if x origin is positive.
+
+    Makes eastern tiles (e.g., 179.7 to 180) adjacent to western tiles
+    (e.g., -180 to -179.8) by shifting to (-180.3 to -180).
+    Only useful for tiles that are part of a dateline-crossing split.
+    """
+    from osgeo import gdal
+
+    ds = gdal.Open(filepath, gdal.GA_Update)
+    gt = list(ds.GetGeoTransform())
+    if gt[0] > 0:
+        logger.info("Shifting tile {} x origin from {} to {}".format(
+            filepath, gt[0], gt[0] - 360.0
+        ))
+        gt[0] -= 360.0
+        ds.SetGeoTransform(gt)
+    ds.FlushCache()
+    ds = None
+
+
+def _download_single_bbox(
+    output_name,
+    bbox,
+    vrt_filename,
+    keep_egm,
+    xrate,
+    yrate,
+    output_format,
+    output_type,
+):
+    """Download a single bbox from the COP DEM."""
+    from osgeo import gdal
+
     if keep_egm:
         t_srs = s_srs = None
     else:
@@ -55,7 +152,6 @@ def download_and_stitch(
     yres = DEFAULT_RES / yrate
     resamp = "bilinear" if (xrate > 1 or yrate > 1) else "nearest"
 
-    # access_mode = "overwrite" if overwrite else None
     option_dict = dict(
         format=output_format,
         outputBounds=utils.align_bounds_to_pixel_grid(bbox),
@@ -69,29 +165,24 @@ def download_and_stitch(
         warpMemoryLimit=5000,
         warpOptions=["NUM_THREADS=4"],
     )
-    # When converting from geoid to ellipsoid, preserve ocean (value=0) as nodata
-    # COP DEM has ocean areas as 0 (sea level relative to geoid), which would
-    # otherwise become ~geoid_undulation after the vertical datum conversion
+    # Preserve ocean (value=0) as nodata during geoid-to-ellipsoid conversion
     if not keep_egm:
         option_dict["srcNodata"] = 0
         option_dict["dstNodata"] = 0
 
-    # Used the __RETURN_OPTION_LIST__ to get the list of options for debugging
     logger.info("Creating {}".format(output_name))
     logger.info("Fetching remote tiles...")
     try:
         cmd = _gdal_cmd_from_options(vrt_filename, output_name, option_dict)
         logger.info("Running GDAL command:")
         logger.info(cmd)
     except Exception:
-        # Can't form the cli version due to `deepcopy` Pickle error, just skip
         logger.info("Running gdal.Warp with options:")
         logger.info(option_dict)
         pass
-    # Now convert to something GDAL can actually use
+
     option_dict["callback"] = gdal.TermProgress
     gdal.Warp(output_name, vrt_filename, options=gdal.WarpOptions(**option_dict))
-    return
 
 
 def _gdal_cmd_from_options(src, dst, option_dict):

sardem/dem.py

@@ -414,24 +414,124 @@ def main(
             output_format,
         )
 
-    tile_names = list(Tile(*bbox).srtm1_tile_names())
+    # Check for dateline crossing
+    bboxes = utils.check_dateline(bbox)
 
-    d = Downloader(tile_names, data_source=data_source, cache_dir=cache_dir)
-    local_filenames = d.download_all()
+    if len(bboxes) == 1:
+        # No dateline crossing, proceed normally
+        tile_names = list(Tile(*bbox).srtm1_tile_names())
 
-    s = Stitcher(tile_names, filenames=local_filenames, data_source=data_source)
-    stitched_dem = s.load_and_stitch()
+        d = Downloader(tile_names, data_source=data_source, cache_dir=cache_dir)
+        local_filenames = d.download_all()
 
-    # Now create corresponding rsc file for all the tiles
-    rsc_dict_tiles = s.create_dem_rsc()
+        s = Stitcher(tile_names, filenames=local_filenames, data_source=data_source)
+        stitched_dem = s.load_and_stitch()
 
-    logger.info("Cropping stitched DEM to boundaries")
-    stitched_dem = upsample.resample(stitched_dem, rsc_dict_tiles, bbox)
-    rsc_dict = rsc_dict_tiles.copy()
-    rsc_dict["X_FIRST"] = bbox[0]
-    rsc_dict["Y_FIRST"] = bbox[3]
-    rsc_dict["FILE_LENGTH"] = stitched_dem.shape[0]
-    rsc_dict["WIDTH"] = stitched_dem.shape[1]
+        rsc_dict_tiles = s.create_dem_rsc()
+
+        logger.info("Cropping stitched DEM to boundaries")
+        stitched_dem = upsample.resample(stitched_dem, rsc_dict_tiles, bbox)
+        rsc_dict = rsc_dict_tiles.copy()
+        rsc_dict["X_FIRST"] = bbox[0]
+        rsc_dict["Y_FIRST"] = bbox[3]
+        rsc_dict["FILE_LENGTH"] = stitched_dem.shape[0]
+        rsc_dict["WIDTH"] = stitched_dem.shape[1]
+    else:
+        # Dateline crossing: download each sub-bbox, merge via GDAL VRT
+        import tempfile
+
+        logger.info(
+            "Dateline crossing detected, downloading {} separate regions".format(
+                len(bboxes)
+            )
+        )
+        utils._gdal_installed_correctly()
+        from osgeo import gdal, osr
+
+        temp_files = []
+        with tempfile.TemporaryDirectory() as tmpdir:
+            for idx, sub_bbox in enumerate(bboxes):
+                logger.info(
+                    "Processing region {} of {}".format(idx + 1, len(bboxes))
+                )
+
+                tile_names = list(Tile(*sub_bbox).srtm1_tile_names())
+                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
+                )
+                dem_part = s.load_and_stitch()
+
+                rsc_dict_part = s.create_dem_rsc()
+                dem_part = upsample.resample(dem_part, rsc_dict_part, sub_bbox)
+
+                # Write as GeoTIFF for VRT merging
+                temp_tif = os.path.join(tmpdir, "dem_part_{}.tif".format(idx))
+                temp_files.append(temp_tif)
+
+                x_first = sub_bbox[0]
+                y_first = sub_bbox[3]
+                x_step = (sub_bbox[2] - sub_bbox[0]) / dem_part.shape[1]
+                y_step = -(sub_bbox[3] - sub_bbox[1]) / dem_part.shape[0]
+
+                # Shift western tile to 0-360 range so GDAL sees tiles as adjacent
+                if x_first < 0:
+                    x_first += 360.0
+
+                driver = gdal.GetDriverByName("GTiff")
+                ds = driver.Create(
+                    temp_tif,
+                    dem_part.shape[1],
+                    dem_part.shape[0],
+                    1,
+                    gdal.GDT_Int16
+                    if data_source != "NASA_WATER"
+                    else gdal.GDT_Byte,
+                )
+                ds.SetGeoTransform([x_first, x_step, 0, y_first, 0, y_step])
+
+                srs = osr.SpatialReference()
+                srs.ImportFromEPSG(4326)
+                ds.SetProjection(srs.ExportToWkt())
+
+                ds.GetRasterBand(1).WriteArray(dem_part)
+                ds.FlushCache()
+                ds = None
+
+            # Merge sub-regions via VRT
+            logger.info("Merging {} regions into VRT".format(len(temp_files)))
+            vrt_temp = os.path.join(tmpdir, "merged.vrt")
+            gdal.BuildVRT(vrt_temp, temp_files)
+
+            ds = gdal.Open(vrt_temp)
+            stitched_dem = ds.GetRasterBand(1).ReadAsArray()
+            geotransform = ds.GetGeoTransform()
+            ds = None
+
+            rsc_dict = collections.OrderedDict.fromkeys(RSC_KEYS)
+            rsc_dict.update(
+                {
+                    "X_UNIT": "degrees",
+                    "Y_UNIT": "degrees",
+                    "Z_OFFSET": 0,
+                    "Z_SCALE": 1,
+                    "PROJECTION": "LL",
+                }
+            )
+            x_first = geotransform[0]
+            # Wrap back to -180..180 if needed
+            if x_first > 180:
+                x_first -= 360
+            rsc_dict["X_FIRST"] = x_first
+            rsc_dict["Y_FIRST"] = geotransform[3]
+            rsc_dict["X_STEP"] = geotransform[1]
+            rsc_dict["Y_STEP"] = geotransform[5]
+            rsc_dict["FILE_LENGTH"] = stitched_dem.shape[0]
+            rsc_dict["WIDTH"] = stitched_dem.shape[1]
 
     rsc_filename = output_name + ".rsc"
 

sardem/tests/test_utils.py

@@ -9,3 +9,66 @@ def test_shift_integer_bbox():
     hp = DEFAULT_RES / 2
     expected = [-156.0 - hp, 19.0 + hp, -155.0 - hp, 20.0 + hp]
     assert utils.shift_integer_bbox(bbox) == pytest.approx(expected)
+
+
+class TestCheckDateline:
+    """Tests for the check_dateline function."""
+
+    def test_no_dateline_crossing(self):
+        bbox = (-156.0, 19.0, -155.0, 20.0)
+        result = utils.check_dateline(bbox)
+        assert len(result) == 1
+        assert result[0] == pytest.approx(bbox)
+
+    def test_no_dateline_crossing_positive_lon(self):
+        bbox = (100.0, -10.0, 110.0, 10.0)
+        result = utils.check_dateline(bbox)
+        assert len(result) == 1
+        assert result[0] == pytest.approx(bbox)
+
+    def test_dateline_crossing_standard(self):
+        # 170E to 170W = from 170 to -170 (crossing 180)
+        bbox = (170.0, -10.0, -170.0, 10.0)
+        result = utils.check_dateline(bbox)
+        assert len(result) == 2
+
+        result = sorted(result, key=lambda x: x[0])
+
+        # Western part: -180 to -170
+        assert result[0][0] == pytest.approx(-180.0)
+        assert result[0][2] == pytest.approx(-170.0)
+        assert result[0][1] == pytest.approx(-10.0)
+        assert result[0][3] == pytest.approx(10.0)
+
+        # Eastern part: 170 to 180
+        assert result[1][0] == pytest.approx(170.0)
+        assert result[1][2] == pytest.approx(180.0)
+        assert result[1][1] == pytest.approx(-10.0)
+        assert result[1][3] == pytest.approx(10.0)
+
+    def test_dateline_crossing_wide_span(self):
+        # 160E to 160W = 40 degrees, crossing dateline
+        bbox = (160.0, -5.0, -160.0, 5.0)
+        result = utils.check_dateline(bbox)
+        assert len(result) == 2
+
+        result = sorted(result, key=lambda x: x[0])
+
+        assert result[0][0] == pytest.approx(-180.0)
+        assert result[0][2] == pytest.approx(-160.0)
+
+        assert result[1][0] == pytest.approx(160.0)
+        assert result[1][2] == pytest.approx(180.0)
+
+    def test_bbox_at_dateline_not_crossing(self):
+        # 175E to 180 (right at dateline but not crossing)
+        bbox = (175.0, -10.0, 180.0, 10.0)
+        result = utils.check_dateline(bbox)
+        assert len(result) == 1
+        assert result[0] == pytest.approx(bbox)
+
+    def test_bbox_west_of_dateline_not_crossing(self):
+        bbox = (-180.0, -10.0, -170.0, 10.0)
+        result = utils.check_dateline(bbox)
+        assert len(result) == 1
+        assert result[0] == pytest.approx(bbox)