Added additional reprojection work for shakemaps that cross the dateline

shakyground2/shakemap.py

@@ -359,6 +359,7 @@ class Shakemap(object):
             )
             * rasterio.transform.Affine.scale(spcx, spcy)
         )
+
         # Export to file
         kwargs0 = {
             "height": results.shape[0],
@@ -366,9 +367,41 @@ class Shakemap(object):
             "count": 1,
             "dtype": results.dtype,
             "transform": transform,
+            "src_transform": transform,
+            "dst_transform": transform,
             "crs": "+proj=latlong",
             "compress": "lzw",
         }
+        if self.site_model.cross_antimeridian:
+            # The code so far defines a transformation on the left lower
+            # corner of the dataset - with positive spacings for each
+            # cell.
+            #
+            # However, if we want to display those cross antimerdian shakemaps
+            # in a geoserver, it causes some trouble on displaying them.
+            # It is possible to run gdal warp on this, but we would favor a
+            # solution that doesn't need another (external) processing step.
+            #
+            # The stuff that the gdal warp fixes in that the transformation
+            # is then defined in the upper left corner. This also changes
+            # the spacing, as the spacing in y direction now needs to be
+            # negative.
+            geoserver_transform = (
+                rasterio.transform.Affine.translation(
+                    self.site_model.bbox_properties["bbox"][0] - (spcx / 2.0),
+                    self.site_model.bbox_properties["bbox"][3] + (spcy / 2.0),
+                )
+                * rasterio.transform.Affine.scale(spcx, spcy * -1.0)
+            )
+            # We need to create those tiffs with the changed transformation
+            kwargs0["transform"] = geoserver_transform
+
+            # However, the data itself still follows the old ordering
+            # (starting from the left lower corner).
+            # As we will run a reproection then anyway, the ordering can
+            # be fixed there as well.
+            kwargs0["dst_transform"] = geoserver_transform
+
         if filename:
             with rasterio.open(filename, "w", "GTiff", **kwargs0) as dst:
                 if self.site_model.cross_antimeridian:
@@ -376,9 +409,9 @@ class Shakemap(object):
                     reproject(
                         results,
                         reproj_results,
-                        src_transform=kwargs0["transform"],
+                        src_transform=kwargs0["src_transform"],
                         src_crs=kwargs0["crs"],
-                        dst_transform=kwargs0["transform"],
+                        dst_transform=kwargs0["dst_transform"],
                         dst_crs=kwargs0["crs"],
                         kwargs={"CENTER_LONG": 180},
                     )
@@ -394,9 +427,9 @@ class Shakemap(object):
                     reproject(
                         results,
                         reproj_results,
-                        src_transform=kwargs0["transform"],
+                        src_transform=kwargs0["src_transform"],
                         src_crs=kwargs0["crs"],
-                        dst_transform=kwargs0["transform"],
+                        dst_transform=kwargs0["dst_transform"],
                         dst_crs=kwargs0["crs"],
                         kwargs={"CENTER_LONG": 180},
                     )