From 29a8fbd8d52b69be41ebef20a4da828ecd1c01bf Mon Sep 17 00:00:00 2001 From: Nils Brinckmann Date: Mon, 1 Nov 2021 16:32:01 +0100 Subject: [PATCH 1/3] Use orientation from upper left corner for all geotiffs --- shakyground2/shakemap.py | 53 ++++++++++++---------------------------- 1 file changed, 15 insertions(+), 38 deletions(-) diff --git a/shakyground2/shakemap.py b/shakyground2/shakemap.py index 9f999fa..b508cf1 100644 --- a/shakyground2/shakemap.py +++ b/shakyground2/shakemap.py @@ -354,12 +354,21 @@ 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) + # In the esri ascii grid we specify the center of lower left corner as reference + # which is bbox[0], bbox[1] + # For the Geotiff we need to specify the upper left corner of the overall raster + # as reference. + # This is still bbox[0] - spcx/2.0 for the x coordinate, + # but the y coordinate we use the bbox[3] - and add then also the half spacing + # to be at the end of the overall cell (upper end). + # As we know need to decrease the y value with every row that we write, we + # need to use a negative spcy. 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 +378,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 +391,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,9 +398,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}, ) @@ -439,9 +416,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}, ) -- GitLab From e7cc05ee94b16db01ab03dfd1485277d29c8c640 Mon Sep 17 00:00:00 2001 From: Nils Brinckmann Date: Tue, 2 Nov 2021 11:38:37 +0100 Subject: [PATCH 2/3] Tried to improve comments --- shakyground2/shakemap.py | 41 +++++++++++++++++++++++++++++++--------- 1 file changed, 32 insertions(+), 9 deletions(-) diff --git a/shakyground2/shakemap.py b/shakyground2/shakemap.py index b508cf1..5a204c1 100644 --- a/shakyground2/shakemap.py +++ b/shakyground2/shakemap.py @@ -354,15 +354,36 @@ 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) - # In the esri ascii grid we specify the center of lower left corner as reference - # which is bbox[0], bbox[1] - # For the Geotiff we need to specify the upper left corner of the overall raster - # as reference. - # This is still bbox[0] - spcx/2.0 for the x coordinate, - # but the y coordinate we use the bbox[3] - and add then also the half spacing - # to be at the end of the overall cell (upper end). - # As we know need to decrease the y value with every row that we write, we - # need to use a negative spcy. + # 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), @@ -402,6 +423,8 @@ class Shakemap(object): src_crs=kwargs0["crs"], 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) -- GitLab From 889167da6ccf63e97c965d26417225cea1b6910a Mon Sep 17 00:00:00 2001 From: Nils Brinckmann Date: Tue, 2 Nov 2021 11:53:14 +0100 Subject: [PATCH 3/3] Removed trailing whitespace --- shakyground2/shakemap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/shakyground2/shakemap.py b/shakyground2/shakemap.py index 5a204c1..9027080 100644 --- a/shakyground2/shakemap.py +++ b/shakyground2/shakemap.py @@ -368,7 +368,7 @@ class Shakemap(object): # 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. # -- GitLab