Use orientation from upper left corner for all geotiffs

shakyground2/shakemap.py

@@ -354,12 +354,42 @@ class Shakemap(object):
         spcx = self.site_model.bbox_properties["spcx"]
         spcy = self.site_model.bbox_properties["spcy"]
         # Build the transformation object (holds the geospatial information about the raster)
+        # The current grid uses the following orientation:
+        #    results[ 0,  0] => upper left corner  (NW)
+        #    results[ 0, -1] => upper right corner (NE)
+        #    results[-1,  0] => lower left corner  (SW)
+        #    results[-1, -1] => lower right corner (SE)
+        #
+        # So with every row that we have we go into the south and with every value in
+        # the row we go into the east.
+        #
+        # Our bounding box has the following orientation:
+        #    bbox[0] => min x
+        #    bbox[1] => min y
+        #    bbox[2] => max x
+        #    bbox[3] => max y
+        #
+        # Geotiffs (unlike esri ascii rasters) use the upper left corner for the default
+        # orientation.
+        #
+        # So we need:
+        #   bbbox[0], bbox[3]
+        #
+        # When we have those bbox corner as middle of the pixel (as it is for our
+        # shakemap rasters here), we need to add/substract the spacing to be
+        # on the upper left corner of the most upper left pixel).
+        # For the bbox[0] we need to substract half of the x spacing.
+        # For the bbox[3] we need to add half of the y spacing.
+        #
+        # Also note that, due to beginning on the upper left corner, we need to
+        # specify a negative y spacing to the transformation.
+        # (The y values for the next rows after the upper left corner must decrease).
         transform = (
             rasterio.transform.Affine.translation(
                 self.site_model.bbox_properties["bbox"][0] - (spcx / 2.0),
-                self.site_model.bbox_properties["bbox"][1] - (spcy / 2.0),
+                self.site_model.bbox_properties["bbox"][3] + (spcy / 2.0),
             )
-            * rasterio.transform.Affine.scale(spcx, spcy)
+            * rasterio.transform.Affine.scale(spcx, -1 * spcy)
         )
 
         # Export to file
@@ -369,8 +399,6 @@ class Shakemap(object):
             "count": 1,
             "dtype": results.dtype,
             "transform": transform,
-            "src_transform": transform,
-            "dst_transform": transform,
             "crs": "+proj=latlong",
             # As we create a lot of those geotiff files in the earthquake event explorer
             # and want to store them on a map server, we want to make sure that the resulting
@@ -384,36 +412,6 @@ class Shakemap(object):
             "zlevel": 9,
             "predictor": 3,
         }
-        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:
@@ -421,10 +419,12 @@ class Shakemap(object):
                     reproject(
                         results,
                         reproj_results,
-                        src_transform=kwargs0["src_transform"],
+                        src_transform=kwargs0["transform"],
                         src_crs=kwargs0["crs"],
-                        dst_transform=kwargs0["dst_transform"],
+                        dst_transform=kwargs0["transform"],
                         dst_crs=kwargs0["crs"],
+                        # We use the very same transformation, but we need the center
+                        # on 180, in order to not run into dateline issues.
                         kwargs={"CENTER_LONG": 180},
                     )
                     dst.write(reproj_results, 1)
@@ -439,9 +439,9 @@ class Shakemap(object):
                     reproject(
                         results,
                         reproj_results,
-                        src_transform=kwargs0["src_transform"],
+                        src_transform=kwargs0["transform"],
                         src_crs=kwargs0["crs"],
-                        dst_transform=kwargs0["dst_transform"],
+                        dst_transform=kwargs0["transform"],
                         dst_crs=kwargs0["crs"],
                         kwargs={"CENTER_LONG": 180},
                     )