diff --git a/shakyground2/shakemap.py b/shakyground2/shakemap.py index 9f999fa812fd9dcea087fa147173de3dc1b7a649..90270809ba51f1fb66aa511dd2aca2515bea6870 100644 --- a/shakyground2/shakemap.py +++ b/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}, )