From cb349511401e5fc2a59c3e527ea237a0ea79ca10 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 12:49:51 -0700 Subject: [PATCH 01/16] Abstracts away common wps dataset operations for gdal --- api/app/utils/wps_dataset.py | 188 +++++++++++++++++++++++++++++++++++ 1 file changed, 188 insertions(+) create mode 100644 api/app/utils/wps_dataset.py diff --git a/api/app/utils/wps_dataset.py b/api/app/utils/wps_dataset.py new file mode 100644 index 000000000..13262ca71 --- /dev/null +++ b/api/app/utils/wps_dataset.py @@ -0,0 +1,188 @@ +from osgeo import gdal, osr +import numpy as np + +from app.utils.geospatial import GDALResamplingMethod + + +class WPSDataset: + """ + A wrapper around gdal datasets for common operations + """ + + def __init__(self, ds_path: str, band=1, chunk_size=256, size=gdal.GDT_Byte): + self.ds = None + self.ds_path = ds_path + self.band = band + self.chunk_size = chunk_size + self.size = size + + def __enter__(self): + self.ds: gdal.Dataset = gdal.Open(self.ds_path) + return self + + def __exit__(self, *_): + self.ds = None + + def __mul__(self, other): + """ + Multiplies this WPSDataset with the other WPSDataset + + :param other: WPSDataset + :raises ValueError: Raised if this and other WPSDataset have mismatched raster dimensions + :return: a new WPSDataset + """ + # Get raster dimensions + x_size = self.ds.RasterXSize + y_size = self.ds.RasterYSize + + # Check if the dimensions of both rasters match + if x_size != other.ds.RasterXSize or y_size != other.ds.RasterYSize: + raise ValueError("The dimensions of the two rasters do not match.") + + # Get the geotransform and projection from the first raster + geotransform = self.ds.GetGeoTransform() + projection = self.ds.GetProjection() + + # Create the output raster + driver: gdal.Driver = gdal.GetDriverByName("MEM") + out_ds: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, self.size) + + # Set the geotransform and projection + out_ds.SetGeoTransform(geotransform) + out_ds.SetProjection(projection) + + self_band: gdal.Band = self.ds.GetRasterBand(self.band) + other_band: gdal.Band = other.ds.GetRasterBand(self.band) + + # Process in chunks + for y in range(0, y_size, self.chunk_size): + y_chunk_size = min(self.chunk_size, y_size - y) + + for x in range(0, x_size, self.chunk_size): + x_chunk_size = min(self.chunk_size, x_size - x) + + # Read chunks from both rasters + self_chunk = self_band.ReadAsArray(x, y, x_chunk_size, y_chunk_size) + other_chunk = other_band.ReadAsArray(x, y, x_chunk_size, y_chunk_size) + + other_chunk[other_chunk >= 1] = 1 + other_chunk[other_chunk < 1] = 0 + + # Multiply the chunks + self_chunk *= other_chunk + + # Write the result to the output raster + out_ds.GetRasterBand(self.band).WriteArray(self_chunk, x, y) + self_chunk = None + other_chunk = None + + return WPSDataset(ds=out_ds) + + def warp_to_match(self, other, output_path: str, resample_method: GDALResamplingMethod = GDALResamplingMethod.NEAREST_NEIGHBOUR) -> gdal.Dataset: + """ + Warp the dataset to match the extent, pixel size, and projection of the other dataset. + + :param ds_to_match: the reference dataset raster to match the source against + :param output_path: output path of the resulting raster + :param resample_method: gdal resampling algorithm + :return: warped raster dataset + """ + source_geotransform = other.ds.GetGeoTransform() + x_res = source_geotransform[1] + y_res = -source_geotransform[5] + minx = source_geotransform[0] + maxy = source_geotransform[3] + maxx = minx + source_geotransform[1] * other.ds.RasterXSize + miny = maxy + source_geotransform[5] * other.ds.RasterYSize + extent = [minx, miny, maxx, maxy] + + # Warp to match input option parameters + warped_ds = gdal.Warp(output_path, self.ds, dstSRS=other.ds.GetProjection(), outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=resample_method.value) + return WPSDataset(ds=warped_ds) + + def replace_nodata_with(self, new_no_data_value: int): + """ + Modifies the dataset by replacing the nodata value with the supplied one + + :param new_no_data_value: the new nodata value + """ + band: gdal.Band = self.ds.GetRasterBand(self.band) + nodata_value = band.GetNoDataValue() + array = band.ReadAsArray() + + if nodata_value is not None: + modified_array = np.where(array != nodata_value, array, new_no_data_value) + band.WriteArray(modified_array) + band.SetNoDataValue(new_no_data_value) + band.FlushCache() + + def generate_latitude_array(self): + """ + Transforms this dataset to 4326 to compute the latitude coordinates + + :return: array of latitude coordinates + """ + geotransform = self.ds.GetGeoTransform() + projection = self.ds.GetProjection() + + src_srs = osr.SpatialReference() + src_srs.ImportFromWkt(projection) + + x_size = self.ds.RasterXSize + y_size = self.ds.RasterYSize + + tgt_srs = osr.SpatialReference() + tgt_srs.ImportFromEPSG(4326) + + transform = osr.CoordinateTransformation(src_srs, tgt_srs) + + # empty array to store latitude values + latitudes = np.zeros((y_size, x_size)) + + for y in range(y_size): + for x in range(x_size): + x_coord = geotransform[0] + x * geotransform[1] + y * geotransform[2] + y_coord = geotransform[3] + x * geotransform[4] + y * geotransform[5] + + _, lat, _ = transform.TransformPoint(x_coord, y_coord) + + latitudes[y, x] = lat + + return latitudes + + def export_to_geotiff(self, output_path): + """ + Exports the dataset to a geotiff with the given path + + :param output_path: path to export the geotiff to + """ + driver: gdal.Driver = gdal.GetDriverByName("GTiff") + + geotransform = self.ds.GetGeoTransform() + projection = self.ds.GetProjection() + + band: gdal.Band = self.ds.GetRasterBand(self.band) + nodata_value = band.GetNoDataValue() + array = band.ReadAsArray() + + rows, cols = array.shape + output_dataset: gdal.Dataset = driver.Create(output_path, cols, rows, 1, self.size) + output_dataset.SetGeoTransform(geotransform) + output_dataset.SetProjection(projection) + + output_band: gdal.Band = output_dataset.GetRasterBand(self.band) + output_band.WriteArray(array) + if nodata_value is not None: + output_band.SetNoDataValue(nodata_value) + + output_band.FlushCache() + output_dataset = None + del output_dataset + output_band = None + del output_band + + def as_gdal_ds(self) -> gdal.Dataset: + return self.ds + + def close(self): + self.ds = None From ce446468ce0c96066ef7d3caeee8a5e366fdf859 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 15:39:45 -0700 Subject: [PATCH 02/16] size -> datatype --- api/app/utils/wps_dataset.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/api/app/utils/wps_dataset.py b/api/app/utils/wps_dataset.py index 13262ca71..f337737d0 100644 --- a/api/app/utils/wps_dataset.py +++ b/api/app/utils/wps_dataset.py @@ -9,12 +9,12 @@ class WPSDataset: A wrapper around gdal datasets for common operations """ - def __init__(self, ds_path: str, band=1, chunk_size=256, size=gdal.GDT_Byte): + def __init__(self, ds_path: str, band=1, chunk_size=256, datatype=gdal.GDT_Byte): self.ds = None self.ds_path = ds_path self.band = band self.chunk_size = chunk_size - self.size = size + self.datatype = datatype def __enter__(self): self.ds: gdal.Dataset = gdal.Open(self.ds_path) @@ -45,7 +45,7 @@ def __mul__(self, other): # Create the output raster driver: gdal.Driver = gdal.GetDriverByName("MEM") - out_ds: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, self.size) + out_ds: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, self.datatype) # Set the geotransform and projection out_ds.SetGeoTransform(geotransform) @@ -166,7 +166,7 @@ def export_to_geotiff(self, output_path): array = band.ReadAsArray() rows, cols = array.shape - output_dataset: gdal.Dataset = driver.Create(output_path, cols, rows, 1, self.size) + output_dataset: gdal.Dataset = driver.Create(output_path, cols, rows, 1, self.datatype) output_dataset.SetGeoTransform(geotransform) output_dataset.SetProjection(projection) From dd0e4742012e70ffc0ca2c9e369bb48753349b0b Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 16:48:42 -0700 Subject: [PATCH 03/16] Tests --- .vscode/settings.json | 1 + api/app/tests/utils/test_wps_dataset.py | 79 ++++++++++++++++++++++++ api/app/tests/utils/zero_layer.tif | Bin 0 -> 2131094 bytes api/app/utils/wps_dataset.py | 30 ++++++--- 4 files changed, 103 insertions(+), 7 deletions(-) create mode 100644 api/app/tests/utils/test_wps_dataset.py create mode 100644 api/app/tests/utils/zero_layer.tif diff --git a/.vscode/settings.json b/.vscode/settings.json index 6174d53ea..6851d849e 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -69,6 +69,7 @@ "cutline", "CWFIS", "determinates", + "dtype", "excinfo", "fastapi", "FBAN", diff --git a/api/app/tests/utils/test_wps_dataset.py b/api/app/tests/utils/test_wps_dataset.py new file mode 100644 index 000000000..02302b85e --- /dev/null +++ b/api/app/tests/utils/test_wps_dataset.py @@ -0,0 +1,79 @@ +import os +import numpy as np +from osgeo import osr, gdal + +from app.utils.wps_dataset import WPSDataset + +hfi_tif = os.path.join(os.path.dirname(__file__), "snow_masked_hfi20240810.tif") +zero_tif = os.path.join(os.path.dirname(__file__), "zero_layer.tif") + + +def create_test_dataset(filename, width, height, extent, projection, data_type=gdal.GDT_Float32) -> gdal.Dataset: + """ + Create a test GDAL dataset. + """ + # Create a new GDAL dataset + driver: gdal.Driver = gdal.GetDriverByName("MEM") + dataset: gdal.Dataset = driver.Create(filename, width, height, 1, data_type) + + # Set the geotransform + xmin, xmax, ymin, ymax = extent + xres = (xmax - xmin) / width + yres = (ymax - ymin) / height + geotransform = (xmin, xres, 0, ymax, 0, -yres) # Top-left corner + dataset.SetGeoTransform(geotransform) + + # Set the projection + srs = osr.SpatialReference() + srs.ImportFromEPSG(projection) + dataset.SetProjection(srs.ExportToWkt()) + + # Create some test data (e.g., random values) + data = np.random.rand(height, width).astype(np.float32) # Generate random data + + # Write data to the dataset + dataset.GetRasterBand(1).WriteArray(data) + return dataset + + +def test_raster_with_context(): + """ + with opens the dataset and closes after the context ends + """ + with WPSDataset(hfi_tif) as wps_ds: + assert wps_ds.as_gdal_ds() is not None + + assert wps_ds.as_gdal_ds() is None + + +def test_raster_set_no_data_value(): + with WPSDataset(hfi_tif) as wps_ds: + assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == 0 + + wps_ds.replace_nodata_with(-1) + assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == -1 + + +def test_raster_mul(): + with WPSDataset(hfi_tif) as wps_ds, WPSDataset(zero_tif) as zero_ds: + output_ds = wps_ds * zero_ds + raw_ds = output_ds.as_gdal_ds() + output_values = raw_ds.GetRasterBand(1).ReadAsArray() + assert np.all(output_values == 0) + + +def test_raster_warp(): + # Dataset 1: 100x100 pixels, extent in EPSG:4326 + extent1 = (-10, 10, -10, 10) # xmin, xmax, ymin, ymax + wgs_84_ds = create_test_dataset("test_dataset_1.tif", 100, 100, extent1, 4326) + + # Dataset 2: 200x200 pixels, extent in EPSG:3857 + extent2 = (-20037508.34, 20037508.34, -20037508.34, 20037508.34) + mercator_ds = create_test_dataset("test_dataset_2.tif", 200, 200, extent2, 3857) + + with WPSDataset(ds_path=None, ds=wgs_84_ds) as wps1_ds, WPSDataset(ds_path=None, ds=mercator_ds) as wps2_ds: + output_ds: WPSDataset = wps1_ds.warp_to_match(wps2_ds, "/vsimem/test.tif") + assert output_ds.as_gdal_ds().GetProjection() == wps2_ds.as_gdal_ds().GetProjection() + assert output_ds.as_gdal_ds().GetGeoTransform() == wps2_ds.as_gdal_ds().GetGeoTransform() + assert output_ds.as_gdal_ds().RasterXSize == wps2_ds.as_gdal_ds().RasterXSize + assert output_ds.as_gdal_ds().RasterYSize == wps2_ds.as_gdal_ds().RasterYSize diff --git a/api/app/tests/utils/zero_layer.tif b/api/app/tests/utils/zero_layer.tif new file mode 100644 index 0000000000000000000000000000000000000000..3e20511b4cb1aa849ef220ef16b2ba02b49edd4c GIT binary patch literal 2131094 zcmeI!Yiwk9c?R%j#`gGr-SyaBUz3O-X$WlB$*`aRYh`0tNVb7pObUW%18bmUld=RW z5!)stbX*XZwnRoC%0)y&MHI4&s6l`NL12ME1%a?Af)YV%$}LR=G7Sjm%&du)4=tkV zm#V)x(s`frKkxs{_$STyJ#^>~^3W}F8BPjqDJ1D>Ys8&=PcbR(JUE1i>v z3+-8IEY=!(N6%R=@8rf~GSwiXQflPwM!s*g zEmt*ir8_yj+Ly*s?r-EXtHsPV@~;~CwRaaY-pDgOjguDv0t5&UAn^ZN;J~{Nv;cln ze}~3{=RNj;X3Nlh^9#ReEwm=1k)*vb8#RxeT5LbF+MGF>14$wG4i&o2td<&6NU3%G z<_o&JPOlcynch}yOf!AWzF*%~+^KD?+LiWJz375cU$eiM&Q@(pU#k|<->QxIs8!oD z*s2|8Ki=)nIM50QuKQH0K6-uR?0ReQ?CGILet6gP-6xyt2kJ+!UAXGliR+J_ym{`5 zn{Hgb>EwxHH#F1BUwrED!omIXr;Kz7|Umpfk9mr^7DYa`#c+Lo&tdHjyK zoL=oqS0kTp+4gKdu|12|v}f&Q z?b*7$J=M3jr*@`2o1bq_L{f;Q%bE`%I4uxDlaW%<+f5b-d0NabScZ9D`ow`Qc91MvNYV8`bC}D zUg*r?^EaMIG?@H;!tT@#bZ7J0?o?jcot4*jXXCxyDSxax%U|iv`VYEO zYU|0;WKZhzJ=uO%PZnR&leJSl*}AJI)q8tV`)p4(ztxks|e+exyIOFZE~hd;O_oAS>eo z+1NXf@}Yq&9~;Q}tph3j%|MnuFp&DE2eSQ*fh_)HAZxvY*}7mb)yoG{J35%n8wOMP zi@~h??O-K7|in4VAg*%n9|r#mYy(_`qPK9ecezNZyCzkUkzpJ?}k!+U?{bJ z9Lnau52eyGoR#y3v+?BNl#dK&`4z)izjHXHcMfOi!^5e6aX8!mJe>SP7>}a;0IGXA;qp7`YG@G}Nrt zk7em|W2rwlmhJx-%i{2O)-D>)*1~wI&mT|ir)dcT{e-Wt0z)FK9TJ|o5yW z?#c4&_hkLOds6z?o-BQ3PwGF|lkMVU7AGgOHb0rIXHBO1lF8IgO=k11$yDy0%*tmc zv+=FTl(#3dJTR5@%2Y~Eo66Gjrc%FYD%-15S$xk_);>CwtuIfd`u(ZY3e(w~m`-Ki zbXJ}@osH|KQ+~~Kmft#^^?RmM`pk5ezB!%xPp7lpHi)Nh^5_FvCt@dLA2`}AzKzA>BX zPi9l=oy+D0bE#ZDmzAS)*|_2C9FPO=KG2*89<8m5|N8Xcd5?Xdd3Wf(`Gwzf66jEwkzxjghuG6c9bf&je8@tcpm(cjwphkGCW+b(~?Z|mki0UylgkpKVy literal 0 HcmV?d00001 diff --git a/api/app/utils/wps_dataset.py b/api/app/utils/wps_dataset.py index f337737d0..382884482 100644 --- a/api/app/utils/wps_dataset.py +++ b/api/app/utils/wps_dataset.py @@ -1,3 +1,4 @@ +from typing import Optional from osgeo import gdal, osr import numpy as np @@ -9,15 +10,16 @@ class WPSDataset: A wrapper around gdal datasets for common operations """ - def __init__(self, ds_path: str, band=1, chunk_size=256, datatype=gdal.GDT_Byte): - self.ds = None + def __init__(self, ds_path: Optional[str], ds=None, band=1, chunk_size=256, datatype=gdal.GDT_Byte): + self.ds = ds self.ds_path = ds_path self.band = band self.chunk_size = chunk_size self.datatype = datatype def __enter__(self): - self.ds: gdal.Dataset = gdal.Open(self.ds_path) + if self.ds is None: + self.ds: gdal.Dataset = gdal.Open(self.ds_path) return self def __exit__(self, *_): @@ -64,6 +66,10 @@ def __mul__(self, other): # Read chunks from both rasters self_chunk = self_band.ReadAsArray(x, y, x_chunk_size, y_chunk_size) other_chunk = other_band.ReadAsArray(x, y, x_chunk_size, y_chunk_size) + wider_type = np.promote_types(self_chunk.dtype, other_chunk.dtype) + + self_chunk = self_chunk.astype(wider_type) + other_chunk = other_chunk.astype(wider_type) other_chunk[other_chunk >= 1] = 1 other_chunk[other_chunk < 1] = 0 @@ -76,9 +82,9 @@ def __mul__(self, other): self_chunk = None other_chunk = None - return WPSDataset(ds=out_ds) + return WPSDataset(ds_path=None, ds=out_ds) - def warp_to_match(self, other, output_path: str, resample_method: GDALResamplingMethod = GDALResamplingMethod.NEAREST_NEIGHBOUR) -> gdal.Dataset: + def warp_to_match(self, other, output_path: str, resample_method: GDALResamplingMethod = GDALResamplingMethod.NEAREST_NEIGHBOUR): """ Warp the dataset to match the extent, pixel size, and projection of the other dataset. @@ -97,8 +103,18 @@ def warp_to_match(self, other, output_path: str, resample_method: GDALResampling extent = [minx, miny, maxx, maxy] # Warp to match input option parameters - warped_ds = gdal.Warp(output_path, self.ds, dstSRS=other.ds.GetProjection(), outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=resample_method.value) - return WPSDataset(ds=warped_ds) + warped_ds = gdal.Warp( + output_path, + self.ds, + options=gdal.WarpOptions( + dstSRS=other.ds.GetProjection(), + outputBounds=extent, + xRes=x_res, + yRes=y_res, + resampleAlg=resample_method.value, + ), + ) + return WPSDataset(ds_path=None, ds=warped_ds) def replace_nodata_with(self, new_no_data_value: int): """ From 9fde9d9fd9ece5d15e1eb0f171fd51b132980fdd Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 16:50:44 -0700 Subject: [PATCH 04/16] Cleanup test datasets --- api/app/tests/utils/test_wps_dataset.py | 3 +++ api/app/utils/wps_dataset.py | 3 --- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/api/app/tests/utils/test_wps_dataset.py b/api/app/tests/utils/test_wps_dataset.py index 02302b85e..f12fd3538 100644 --- a/api/app/tests/utils/test_wps_dataset.py +++ b/api/app/tests/utils/test_wps_dataset.py @@ -77,3 +77,6 @@ def test_raster_warp(): assert output_ds.as_gdal_ds().GetGeoTransform() == wps2_ds.as_gdal_ds().GetGeoTransform() assert output_ds.as_gdal_ds().RasterXSize == wps2_ds.as_gdal_ds().RasterXSize assert output_ds.as_gdal_ds().RasterYSize == wps2_ds.as_gdal_ds().RasterYSize + + wgs_84_ds = None + mercator_ds = None diff --git a/api/app/utils/wps_dataset.py b/api/app/utils/wps_dataset.py index 382884482..14619105f 100644 --- a/api/app/utils/wps_dataset.py +++ b/api/app/utils/wps_dataset.py @@ -199,6 +199,3 @@ def export_to_geotiff(self, output_path): def as_gdal_ds(self) -> gdal.Dataset: return self.ds - - def close(self): - self.ds = None From 5b64343d021c6cb2493fcf4aa842996ef3966185 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 17:02:42 -0700 Subject: [PATCH 05/16] update dataset access mode for replacing no datavalue --- .../tests/utils/snow_masked_hfi20240810.tif | Bin 537771 -> 538715 bytes api/app/tests/utils/test_wps_dataset.py | 2 +- api/app/utils/wps_dataset.py | 7 ++++--- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/api/app/tests/utils/snow_masked_hfi20240810.tif b/api/app/tests/utils/snow_masked_hfi20240810.tif index 7e7189c9351dd87d2d96405a31c285ec1d7b9ec8..46f8091afa5213fbf0b8af3dd07431c409836ac2 100644 GIT binary patch delta 142 zcmZ4eQsMRo1$IwQErvBF9F6R)?2N7KOs(w9t?VqV?5wTqY_068(YTxZHLu~gAT#+-wpsfCAuje&u|1jt_B&dlHgWTybx8`^moIDqUwK=#IVCI)#R k+s1-}p>YWl12d4VYX~HOVB&+8i62;) Date: Wed, 16 Oct 2024 17:18:49 -0700 Subject: [PATCH 06/16] remove problem test --- api/app/tests/utils/test_wps_dataset.py | 30 ++++++++++++++++++------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/api/app/tests/utils/test_wps_dataset.py b/api/app/tests/utils/test_wps_dataset.py index 026540905..aa2afc2c9 100644 --- a/api/app/tests/utils/test_wps_dataset.py +++ b/api/app/tests/utils/test_wps_dataset.py @@ -46,14 +46,6 @@ def test_raster_with_context(): assert wps_ds.as_gdal_ds() is None -def test_raster_set_no_data_value(): - with WPSDataset(hfi_tif, access=gdal.GA_Update) as wps_ds: - assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == 0 - - wps_ds.replace_nodata_with(-1) - assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == -1 - - def test_raster_mul(): with WPSDataset(hfi_tif) as wps_ds, WPSDataset(zero_tif) as zero_ds: output_ds = wps_ds * zero_ds @@ -80,3 +72,25 @@ def test_raster_warp(): wgs_84_ds = None mercator_ds = None + + +def test_generate_latitude_array(): + driver: gdal.Driver = gdal.GetDriverByName("MEM") + dataset: gdal.Dataset = driver.Create("test_lat.tif", 3, 3, 1, gdal.GDT_Float32) + + # Set a simple geotransform (top-left corner at (0, 0), 1x1 pixel size, no rotation) + geotransform = (0, 1, 0, 0, 0, -1) + dataset.SetGeoTransform(geotransform) + + # Set projection to EPSG:3857 (Web Mercator) + srs = osr.SpatialReference() + srs.ImportFromEPSG(3857) # Use a known projection + dataset.SetProjection(srs.ExportToWkt()) + + with WPSDataset(ds_path=None, ds=dataset) as lat_ds: + latitudes = lat_ds.generate_latitude_array() + + # Define the expected latitude values based on the transformation logic + expected_latitudes = np.array([[0, 0, 0], [-1, -1, -1], [-2, -2, -2]]) + + np.testing.assert_array_almost_equal(latitudes, expected_latitudes) From c68db6df1025f89fe4a10b954703025b7678e0a5 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 17:28:57 -0700 Subject: [PATCH 07/16] adjust tests --- .vscode/settings.json | 1 + api/app/tests/utils/test_wps_dataset.py | 27 +++---------------------- 2 files changed, 4 insertions(+), 24 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 6851d849e..2c0821110 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -60,6 +60,7 @@ "allclose", "anyio", "APCP", + "astype", "asyncpg", "autoflush", "Behaviour", diff --git a/api/app/tests/utils/test_wps_dataset.py b/api/app/tests/utils/test_wps_dataset.py index aa2afc2c9..135c4d964 100644 --- a/api/app/tests/utils/test_wps_dataset.py +++ b/api/app/tests/utils/test_wps_dataset.py @@ -29,10 +29,11 @@ def create_test_dataset(filename, width, height, extent, projection, data_type=g dataset.SetProjection(srs.ExportToWkt()) # Create some test data (e.g., random values) - data = np.random.rand(height, width).astype(np.float32) # Generate random data + rng = np.random.default_rng(seed=42) # Reproducible random generator + random_data = rng.random((height, width)).astype(np.float32) # Write data to the dataset - dataset.GetRasterBand(1).WriteArray(data) + dataset.GetRasterBand(1).WriteArray(random_data) return dataset @@ -72,25 +73,3 @@ def test_raster_warp(): wgs_84_ds = None mercator_ds = None - - -def test_generate_latitude_array(): - driver: gdal.Driver = gdal.GetDriverByName("MEM") - dataset: gdal.Dataset = driver.Create("test_lat.tif", 3, 3, 1, gdal.GDT_Float32) - - # Set a simple geotransform (top-left corner at (0, 0), 1x1 pixel size, no rotation) - geotransform = (0, 1, 0, 0, 0, -1) - dataset.SetGeoTransform(geotransform) - - # Set projection to EPSG:3857 (Web Mercator) - srs = osr.SpatialReference() - srs.ImportFromEPSG(3857) # Use a known projection - dataset.SetProjection(srs.ExportToWkt()) - - with WPSDataset(ds_path=None, ds=dataset) as lat_ds: - latitudes = lat_ds.generate_latitude_array() - - # Define the expected latitude values based on the transformation logic - expected_latitudes = np.array([[0, 0, 0], [-1, -1, -1], [-2, -2, -2]]) - - np.testing.assert_array_almost_equal(latitudes, expected_latitudes) From d00fb4cd452ada3dcadd262896d7921c0ede7630 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 18:03:17 -0700 Subject: [PATCH 08/16] move to geospatial module --- .vscode/settings.json | 1 + .../{geospatial.py => geospatial/__init__.py} | 8 +- api/app/{utils => geospatial}/wps_dataset.py | 0 .../snow_masked_hfi20240810.tif | Bin .../{utils => geospatial}/test_wps_dataset.py | 43 ++++++- .../{utils => geospatial}/zero_layer.tif | Bin api/app/tests/utils/test_geospatial.py | 108 ------------------ 7 files changed, 44 insertions(+), 116 deletions(-) rename api/app/{geospatial.py => geospatial/__init__.py} (80%) rename api/app/{utils => geospatial}/wps_dataset.py (100%) rename api/app/tests/{utils => geospatial}/snow_masked_hfi20240810.tif (100%) rename api/app/tests/{utils => geospatial}/test_wps_dataset.py (61%) rename api/app/tests/{utils => geospatial}/zero_layer.tif (100%) delete mode 100644 api/app/tests/utils/test_geospatial.py diff --git a/.vscode/settings.json b/.vscode/settings.json index 2c0821110..4355543cc 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -113,6 +113,7 @@ "PRIMEM", "PROJCS", "pydantic", + "pyproj", "RDPS", "reduxjs", "reproject", diff --git a/api/app/geospatial.py b/api/app/geospatial/__init__.py similarity index 80% rename from api/app/geospatial.py rename to api/app/geospatial/__init__.py index 89fee7761..8646ed384 100644 --- a/api/app/geospatial.py +++ b/api/app/geospatial/__init__.py @@ -1,5 +1,5 @@ -""" Geospatial things -""" +"""Module for geospatial logic""" + from typing import Final from pyproj import CRS @@ -8,7 +8,7 @@ # BCGOV standard is to store everything in NAD83 / BC Albers (EPSG:3005). NAD83_BC_ALBERS: Final = 3005 # NAD 83 alone (EPSG:4269), uses geographic coordinates. -NAD83: Final = 'epsg:4269' +NAD83: Final = "epsg:4269" NAD83_CRS: Final = CRS(NAD83) # De facto standard is to expose data in WGS84 (EPSG:4326). -WGS84: Final = 'epsg:4326' +WGS84: Final = "epsg:4326" diff --git a/api/app/utils/wps_dataset.py b/api/app/geospatial/wps_dataset.py similarity index 100% rename from api/app/utils/wps_dataset.py rename to api/app/geospatial/wps_dataset.py diff --git a/api/app/tests/utils/snow_masked_hfi20240810.tif b/api/app/tests/geospatial/snow_masked_hfi20240810.tif similarity index 100% rename from api/app/tests/utils/snow_masked_hfi20240810.tif rename to api/app/tests/geospatial/snow_masked_hfi20240810.tif diff --git a/api/app/tests/utils/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py similarity index 61% rename from api/app/tests/utils/test_wps_dataset.py rename to api/app/tests/geospatial/test_wps_dataset.py index 135c4d964..e9023c432 100644 --- a/api/app/tests/utils/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -1,14 +1,15 @@ import os import numpy as np from osgeo import osr, gdal +import pytest -from app.utils.wps_dataset import WPSDataset +from app.geospatial.wps_dataset import WPSDataset hfi_tif = os.path.join(os.path.dirname(__file__), "snow_masked_hfi20240810.tif") zero_tif = os.path.join(os.path.dirname(__file__), "zero_layer.tif") -def create_test_dataset(filename, width, height, extent, projection, data_type=gdal.GDT_Float32) -> gdal.Dataset: +def create_test_dataset(filename, width, height, extent, projection, data_type=gdal.GDT_Float32, fill_value=None) -> gdal.Dataset: """ Create a test GDAL dataset. """ @@ -32,8 +33,13 @@ def create_test_dataset(filename, width, height, extent, projection, data_type=g rng = np.random.default_rng(seed=42) # Reproducible random generator random_data = rng.random((height, width)).astype(np.float32) - # Write data to the dataset - dataset.GetRasterBand(1).WriteArray(random_data) + if fill_value is None: + # Write data to the dataset + dataset.GetRasterBand(1).WriteArray(random_data) + else: + fill_data = np.full_like(random_data, fill_value) + dataset.GetRasterBand(1).WriteArray(fill_data) + return dataset @@ -55,6 +61,35 @@ def test_raster_mul(): assert np.all(output_values == 0) +def test_raster_mul_identity(): + extent = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax + ds_1 = create_test_dataset("test_dataset_1.tif", 1, 1, extent, 4326, data_type=gdal.GDT_Byte, fill_value=2) + ds_2 = create_test_dataset("test_dataset_2.tif", 1, 1, extent, 4326, data_type=gdal.GDT_Byte, fill_value=1) + + with WPSDataset(ds_path=None, ds=ds_1) as wps1_ds, WPSDataset(ds_path=None, ds=ds_2) as wps2_ds: + output_ds = wps1_ds * wps2_ds + output_values = output_ds.as_gdal_ds().GetRasterBand(1).ReadAsArray() + left_side_values = wps1_ds.as_gdal_ds().GetRasterBand(1).ReadAsArray() + assert np.all(output_values == left_side_values) == True + + +def test_raster_mul_wrong_dimensions(): + # Dataset 1: 100x100 pixels, extent in EPSG:4326 + extent1 = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax + wgs_84_ds = create_test_dataset("test_dataset_1.tif", 1, 1, extent1, 4326) + + # Dataset 2: 200x200 pixels, extent in EPSG:3857 + extent2 = (-2, 2, -2, 2) + mercator_ds = create_test_dataset("test_dataset_2.tif", 2, 2, extent2, 3857) + + with pytest.raises(ValueError): + with WPSDataset(ds_path=None, ds=wgs_84_ds) as wps1_ds, WPSDataset(ds_path=None, ds=mercator_ds) as wps2_ds: + _ = wps1_ds * wps2_ds + + wgs_84_ds = None + mercator_ds = None + + def test_raster_warp(): # Dataset 1: 100x100 pixels, extent in EPSG:4326 extent1 = (-10, 10, -10, 10) # xmin, xmax, ymin, ymax diff --git a/api/app/tests/utils/zero_layer.tif b/api/app/tests/geospatial/zero_layer.tif similarity index 100% rename from api/app/tests/utils/zero_layer.tif rename to api/app/tests/geospatial/zero_layer.tif diff --git a/api/app/tests/utils/test_geospatial.py b/api/app/tests/utils/test_geospatial.py deleted file mode 100644 index 7a88b3f70..000000000 --- a/api/app/tests/utils/test_geospatial.py +++ /dev/null @@ -1,108 +0,0 @@ -import os -import pytest -from osgeo import gdal -import numpy as np - -from app.utils.geospatial import raster_mul, warp_to_match_raster - -fixture_path = os.path.join(os.path.dirname(__file__), "snow_masked_hfi20240810.tif") - - -def get_test_tpi_raster(hfi_ds: gdal.Dataset, fill_value: int): - # Get raster dimensions - x_size = hfi_ds.RasterXSize - y_size = hfi_ds.RasterYSize - - # Get the geotransform and projection from the first raster - geotransform = hfi_ds.GetGeoTransform() - projection = hfi_ds.GetProjection() - - # Create the output raster - driver = gdal.GetDriverByName("MEM") - out_ds: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, gdal.GDT_Byte) - - # Set the geotransform and projection - out_ds.SetGeoTransform(geotransform) - out_ds.SetProjection(projection) - - filler_data = hfi_ds.GetRasterBand(1).ReadAsArray() - tpi_data = np.full_like(filler_data, fill_value) - - # Write the modified data to the new raster - out_band = out_ds.GetRasterBand(1) - out_band.SetNoDataValue(0) - out_band.WriteArray(tpi_data) - return out_ds - - -def get_tpi_raster_wrong_shape(): - driver = gdal.GetDriverByName("MEM") - out_ds: gdal.Dataset = driver.Create("memory", 1, 1, 1, gdal.GDT_Byte) - out_band = out_ds.GetRasterBand(1) - out_band.SetNoDataValue(0) - out_band.WriteArray(np.array([[1]])) - return out_ds - - -def test_zero_case(): - hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) - tpi_ds: gdal.Dataset = get_test_tpi_raster(hfi_ds, 0) - - masked_raster = raster_mul(tpi_ds, hfi_ds) - masked_data = masked_raster.GetRasterBand(1).ReadAsArray() - - assert masked_data.shape == hfi_ds.GetRasterBand(1).ReadAsArray().shape - assert np.all(masked_data == 0) == True - - hfi_ds = None - tpi_ds = None - - -def test_identity_case(): - hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) - tpi_ds: gdal.Dataset = get_test_tpi_raster(hfi_ds, 1) - - masked_raster = raster_mul(tpi_ds, hfi_ds) - masked_data = masked_raster.GetRasterBand(1).ReadAsArray() - hfi_data = hfi_ds.GetRasterBand(1).ReadAsArray() - - # do the simple classification for hfi, pixels >4k are 1 - hfi_data[hfi_data >= 1] = 1 - hfi_data[hfi_data < 1] = 0 - - assert masked_data.shape == hfi_data.shape - assert np.all(masked_data == hfi_data) == True - - hfi_ds = None - tpi_ds = None - - -def test_wrong_dimensions(): - hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) - tpi_ds: gdal.Dataset = get_tpi_raster_wrong_shape() - - with pytest.raises(ValueError): - raster_mul(tpi_ds, hfi_ds) - - hfi_ds = None - tpi_ds = None - - -@pytest.mark.skip(reason="enable once gdal is updated past version 3.4") -def test_warp_to_match_dimension(): - hfi_ds: gdal.Dataset = gdal.Open(fixture_path, gdal.GA_ReadOnly) - tpi_ds: gdal.Dataset = get_tpi_raster_wrong_shape() - - driver = gdal.GetDriverByName("MEM") - out_dataset: gdal.Dataset = driver.Create("memory", hfi_ds.RasterXSize, hfi_ds.RasterYSize, 1, gdal.GDT_Byte) - - warp_to_match_raster(tpi_ds, hfi_ds, out_dataset) - output_data = out_dataset.GetRasterBand(1).ReadAsArray() - hfi_data = hfi_ds.GetRasterBand(1).ReadAsArray() - - assert hfi_data.shape == output_data.shape - assert hfi_ds.RasterXSize == out_dataset.RasterXSize - assert hfi_ds.RasterYSize == out_dataset.RasterYSize - - hfi_ds = None - tpi_ds = None From ac5e04ed9bcadb8213fcacb48d92ac7196cad85e Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 18:09:26 -0700 Subject: [PATCH 09/16] Add back test --- api/app/tests/geospatial/test_wps_dataset.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/api/app/tests/geospatial/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py index e9023c432..8809ba418 100644 --- a/api/app/tests/geospatial/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -40,6 +40,8 @@ def create_test_dataset(filename, width, height, extent, projection, data_type=g fill_data = np.full_like(random_data, fill_value) dataset.GetRasterBand(1).WriteArray(fill_data) + dataset.GetRasterBand(1).SetNoDataValue(0) + return dataset @@ -53,6 +55,16 @@ def test_raster_with_context(): assert wps_ds.as_gdal_ds() is None +def test_raster_set_no_data_value(): + extent = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax + ds_1 = create_test_dataset("test_dataset_no_data_value.tif", 1, 1, extent, 4326, data_type=gdal.GDT_Byte, fill_value=2) + with WPSDataset(ds_path=None, ds=ds_1) as wps_ds: + assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == 0 + + wps_ds.replace_nodata_with(-1) + assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == -1 + + def test_raster_mul(): with WPSDataset(hfi_tif) as wps_ds, WPSDataset(zero_tif) as zero_ds: output_ds = wps_ds * zero_ds From 0fd5c84200f6aa2c09fcb3499a87ac9bf2981e4d Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 18:19:22 -0700 Subject: [PATCH 10/16] add export test --- api/app/tests/geospatial/test_wps_dataset.py | 24 ++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/api/app/tests/geospatial/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py index 8809ba418..1867f184e 100644 --- a/api/app/tests/geospatial/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -2,6 +2,7 @@ import numpy as np from osgeo import osr, gdal import pytest +import tempfile from app.geospatial.wps_dataset import WPSDataset @@ -120,3 +121,26 @@ def test_raster_warp(): wgs_84_ds = None mercator_ds = None + + +def test_export_to_geotiff(): + # Dataset 1: 100x100 pixels, extent in EPSG:4326 + extent1 = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax + ds_1 = create_test_dataset("test_dataset_1.tif", 3, 3, extent1, 4326, data_type=gdal.GDT_Byte, fill_value=1) + + with WPSDataset(ds_path=None, ds=ds_1) as wps_ds: + with tempfile.TemporaryDirectory() as temp_dir: + temp_path = os.path.join(temp_dir, "test_export.tif") + wps_ds.export_to_geotiff(temp_path) + + with WPSDataset(ds_path=temp_path) as exported_ds: + assert wps_ds.as_gdal_ds().GetProjection() == exported_ds.as_gdal_ds().GetProjection() + assert wps_ds.as_gdal_ds().GetGeoTransform() == exported_ds.as_gdal_ds().GetGeoTransform() + assert wps_ds.as_gdal_ds().RasterXSize == exported_ds.as_gdal_ds().RasterXSize + assert wps_ds.as_gdal_ds().RasterYSize == exported_ds.as_gdal_ds().RasterYSize + + original_values = wps_ds.as_gdal_ds().GetRasterBand(1).ReadAsArray() + exported_values = exported_ds.as_gdal_ds().GetRasterBand(1).ReadAsArray() + assert np.all(original_values == exported_values) == True + + ds_1 = None From e8295dd52c10f5489009786846a763e30a069280 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Wed, 16 Oct 2024 18:30:10 -0700 Subject: [PATCH 11/16] Clean up comments --- api/app/tests/geospatial/test_wps_dataset.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/api/app/tests/geospatial/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py index 1867f184e..9fefd9980 100644 --- a/api/app/tests/geospatial/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -87,11 +87,9 @@ def test_raster_mul_identity(): def test_raster_mul_wrong_dimensions(): - # Dataset 1: 100x100 pixels, extent in EPSG:4326 extent1 = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax wgs_84_ds = create_test_dataset("test_dataset_1.tif", 1, 1, extent1, 4326) - # Dataset 2: 200x200 pixels, extent in EPSG:3857 extent2 = (-2, 2, -2, 2) mercator_ds = create_test_dataset("test_dataset_2.tif", 2, 2, extent2, 3857) @@ -124,7 +122,6 @@ def test_raster_warp(): def test_export_to_geotiff(): - # Dataset 1: 100x100 pixels, extent in EPSG:4326 extent1 = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax ds_1 = create_test_dataset("test_dataset_1.tif", 3, 3, extent1, 4326, data_type=gdal.GDT_Byte, fill_value=1) From 0c4e81a4577272d8742e604935c5af4bd89e3d69 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Thu, 17 Oct 2024 11:42:23 -0700 Subject: [PATCH 12/16] Test for same projection and origins --- api/app/geospatial/wps_dataset.py | 8 +++++ api/app/tests/geospatial/test_wps_dataset.py | 33 +++++++++++++++++--- 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/api/app/geospatial/wps_dataset.py b/api/app/geospatial/wps_dataset.py index 215741f1e..62e931704 100644 --- a/api/app/geospatial/wps_dataset.py +++ b/api/app/geospatial/wps_dataset.py @@ -46,6 +46,14 @@ def __mul__(self, other): geotransform = self.ds.GetGeoTransform() projection = self.ds.GetProjection() + # Check if projection matches + if projection != other.ds.GetProjection(): + raise ValueError("The projections of the two rasters do not match.") + + # Check if origin matches + if geotransform[0] != other.ds.GetGeoTransform()[0] or geotransform[3] != other.ds.GetGeoTransform()[3]: + raise ValueError("The origins of the two rasters do not match.") + # Create the output raster driver: gdal.Driver = gdal.GetDriverByName("MEM") out_ds: gdal.Dataset = driver.Create("memory", x_size, y_size, 1, self.datatype) diff --git a/api/app/tests/geospatial/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py index 9fefd9980..56f5af071 100644 --- a/api/app/tests/geospatial/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -87,11 +87,22 @@ def test_raster_mul_identity(): def test_raster_mul_wrong_dimensions(): - extent1 = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax - wgs_84_ds = create_test_dataset("test_dataset_1.tif", 1, 1, extent1, 4326) + extent = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax + wgs_84_ds1 = create_test_dataset("test_dataset_1.tif", 1, 1, extent, 4326) + wgs_84_ds2 = create_test_dataset("test_dataset_2.tif", 2, 2, extent, 4326) + + with pytest.raises(ValueError): + with WPSDataset(ds_path=None, ds=wgs_84_ds1) as wps1_ds, WPSDataset(ds_path=None, ds=wgs_84_ds2) as wps2_ds: + _ = wps1_ds * wps2_ds + + wgs_84_ds1 = None + wgs_84_ds2 = None - extent2 = (-2, 2, -2, 2) - mercator_ds = create_test_dataset("test_dataset_2.tif", 2, 2, extent2, 3857) + +def test_raster_mul_wrong_projections(): + extent = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax + wgs_84_ds = create_test_dataset("test_dataset_1.tif", 1, 1, extent, 4326) + mercator_ds = create_test_dataset("test_dataset_2.tif", 1, 1, extent, 3857) with pytest.raises(ValueError): with WPSDataset(ds_path=None, ds=wgs_84_ds) as wps1_ds, WPSDataset(ds_path=None, ds=mercator_ds) as wps2_ds: @@ -101,6 +112,20 @@ def test_raster_mul_wrong_dimensions(): mercator_ds = None +def test_raster_mul_wrong_origins(): + extent1 = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax + wgs_84_ds1 = create_test_dataset("test_dataset_1.tif", 1, 1, extent1, 4326) + extent2 = (-2, 2, -2, 2) # xmin, xmax, ymin, ymax + wgs_84_ds2 = create_test_dataset("test_dataset_2.tif", 1, 1, extent2, 4326) + + with pytest.raises(ValueError): + with WPSDataset(ds_path=None, ds=wgs_84_ds1) as wps1_ds, WPSDataset(ds_path=None, ds=wgs_84_ds2) as wps2_ds: + _ = wps1_ds * wps2_ds + + wgs_84_ds1 = None + wgs_84_ds2 = None + + def test_raster_warp(): # Dataset 1: 100x100 pixels, extent in EPSG:4326 extent1 = (-10, 10, -10, 10) # xmin, xmax, ymin, ymax From 738c910a5db73eec70f23c4749294cb73d36eed4 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Thu, 17 Oct 2024 16:12:10 -0700 Subject: [PATCH 13/16] test for generate latitude array --- api/app/tests/geospatial/3005_lats.tif | Bin 0 -> 82345 bytes api/app/tests/geospatial/4326_lats.tif | Bin 0 -> 70478 bytes api/app/tests/geospatial/test_wps_dataset.py | 10 ++++++++++ 3 files changed, 10 insertions(+) create mode 100644 api/app/tests/geospatial/3005_lats.tif create mode 100644 api/app/tests/geospatial/4326_lats.tif diff --git a/api/app/tests/geospatial/3005_lats.tif b/api/app/tests/geospatial/3005_lats.tif new file mode 100644 index 0000000000000000000000000000000000000000..1ffec992d4b1d4e88a04b81627c8757e986f421e GIT binary patch literal 82345 zcmeI!F^e2^7{~Ev?~d~zXNW-!LI@jD84HUvK}5Vng@A}C2!b|33P~E1_KFZAXQP6p zr1J&HktzhiS0Gr0KwgD~jR^kt{0B107@i+-_0yQs(CmgBW@{9L*Hz5K4%$5+qWeC+sz#YP+{kNsJGU#gFfcu5=f_h||C~Jfpgi}=GUA`NB7XQh;`DCBou4A!y&G}=e#Do5 zMg0BZBECCUPkM0X(xmb8&57Ih-v9L4+mCMRlm7SS_V#PnFTZze=jT5@`|85(8#|Al zxH@aBu6d)rL%m(w*xx%|{<7HG-`kAiv00Y1i0{wD_UiW6BO3>|ub%l}=biKCmai`_ zTwLy4{o?Yq>$?Y!lxuST_vU{CHgDXAX1);=Xp)lejo=CyahL+KCHV!zTo>;;SOHl` zk1rT30&mg=Ysu?bz#sPRo`4|Bd0LxXc(#}uyJkW?LM3U%jv4TiGj zyp9Q7JF;Rz+jP1VWFaehNr0?gvO0Zv3^WU}kafu4S&)L%A(N_Kj>tk5vdX7mp#LoJ zwt9KKv-D;-%QQn)rV^HfN`|EBNUTvU35@I1M6U&fP9f|#A_UBv4 zP?m2jCKIqDONfOmAr`WPSjZA$AxnsbEFl)MgjmSR#DY|&5|;Ez$kLnPEYl2GnMzpF zDU35@I1M z6AMzAN?6h>Axm$DvrIE&Wh!AwuY@eU8O}1zkd>)~C7}|egjmQDVj)Y2g)AW!vV>U3 z5@I1sh=nX67P5p`$P!{X$zpNB;(TK{k*c%!#$qzTJhFsX$P!{9ONfOmAr`WPSjZA$ zAxnsbtV}FOWh!AwuY@eU8O}1zkd>)~CA|`|^kz8AG(%RV5|;EzPO{jazCXiePNeEA zhLspe$Brz$8O}1zkd>)~CA|`|^kz8AG;@;0#xfhTujE9k&SKw;0W$(knU1Vt@Mne48O_O53z>=0vK_V$p>mx_;D<)wLi7+Y2JA zy$fEtc0pFxf*5Qs=p>6BweM(TU7SeOSw>cnaoax#$ZGF`m#$rq)wLi7+Y2JAy$fEt zc5#x$HoI=qsUWhZv`q_HPNeEA7F`&k>qiY)T?=Bcy&$sMyWpj37bjV4v+FjU3LsFY&DU{{0-Xy~e4WqJ_X++5JfxUd literal 0 HcmV?d00001 diff --git a/api/app/tests/geospatial/4326_lats.tif b/api/app/tests/geospatial/4326_lats.tif new file mode 100644 index 0000000000000000000000000000000000000000..d11d5347dd2fc0974bbf2f9138c0acbbd0700e6f GIT binary patch literal 70478 zcmeIyzf06{90u^ushFjNLKZ}sP*XuPG!%qc4dPIPhBdSVO+`yIMFjo4riO;-4~W(j zA{rc;`U|3=!7+k1p)IiYb6#LUQ|rs)xKCfc_lcJxkbXJ59q^scq*pIW=O-rDo={^K?uT0&c2bE|7_ z?_aoa`pD$wXaA;uek{D7di?$B)ek4n{n*&pI6wRQ#yqp-mk*VF{i5RR{Xh9Zo%VIS zGqTc+%W&`C&kk-MS?=aLqrK18p5DCuV*jgKdEV^niPcar|yORTt{d2LdTj z1X~aRTYx|c6d?k(0D%-JLIi990x3|02-pGyQlJPCumuRDKoKHf3lK~5wHaaq(BiOU<(jPfg(h} z79fxUMTmecKp+K*5CL0&KnfHg0=58w6evOjYyko(P=pBB0t8Z^2obOa2&6y}B47&; zNP!|mz!o5o0!4^`EkGaziVy)?fItcqAp*7lffOi01Z)8UDNuw6*a8Gnpa>DL1qh@- z5h7p<5J-U{M8FmxkOD=BfGt2E1&R;>TYx|c6d?k(0D%-JLIi990x3|02-pGyQlJPC zumuRDKoKHf3lK~5wHaaq(BiOU<(jPfg(h}79fxUMTmecKp+K*5CL0&KnfHg0=58w6evOjYyko( zP=pBB0t8Z^2obOa2&6y}B47&;NP!|mz!o5o0!4^`EkGaziVy)?fItcqAp*7lffOi0 z1Z)8UDNuw6*a8Gnpa>DL1qh@-5h7p<5J-U{M8FmxkOD=BfGt2E1&R;>TYx|c6d?k( P0D%-JLIiBVe-r!yJrQ4s literal 0 HcmV?d00001 diff --git a/api/app/tests/geospatial/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py index 56f5af071..55fc8f134 100644 --- a/api/app/tests/geospatial/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -166,3 +166,13 @@ def test_export_to_geotiff(): assert np.all(original_values == exported_values) == True ds_1 = None + + +def test_latitude_array(): + lats_3005_tif = os.path.join(os.path.dirname(__file__), "3005_lats.tif") + lats_4326_tif = os.path.join(os.path.dirname(__file__), "4326_lats.tif") + with WPSDataset(ds_path=lats_3005_tif) as lats_3005_ds, WPSDataset(ds_path=lats_4326_tif) as lats_4326_ds: + output_ds: WPSDataset = lats_3005_ds.warp_to_match(lats_4326_ds, "/vsimem/test_lats.tif") + original_lats = lats_4326_ds.generate_latitude_array() + warped_lats = output_ds.generate_latitude_array() + assert np.all(original_lats == warped_lats) == True From ad6cd85c06d1af70e5a4ecfb52b69a42cdf05073 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Thu, 17 Oct 2024 17:18:24 -0700 Subject: [PATCH 14/16] Adjust test --- api/app/tests/geospatial/3005_lats.tif | Bin 82345 -> 94001 bytes api/app/tests/geospatial/4326_lats.tif | Bin 70478 -> 80426 bytes api/app/tests/geospatial/test_wps_dataset.py | 4 +++- 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/api/app/tests/geospatial/3005_lats.tif b/api/app/tests/geospatial/3005_lats.tif index 1ffec992d4b1d4e88a04b81627c8757e986f421e..09f69d35eb7f659d0b91cfe1f2b19c2403bf54af 100644 GIT binary patch literal 94001 zcmeI%f2^K$c?R(3aIjimiQ_EzE&oSL}DM->tsHO?dVu_6{ z8mm$o$zU{rt;->5YOEs@cZ5umhBCV}p$p4!MbJuD>Z)L?RltLwX#Ey@OLLX$T-S|h zp6`KIg7sGlzgPTKGfM-{(a36oxeH!dbs3&@|l^whY*5B=vKUi@%iP4mV3n&!<1n&z7`hMGlZ3^jY+G1Odn`A~Ds$k^t0U;NQF z=I`Ev|M%F;H?COvjuRJb{MMyczWwU+7DRLC+0XKin&$ZOX&(CL@I3DL@oP^mrJ)nX zubtMM+)OJU{!nxM>zkQ}=3O(6e|Dbn`q2d!%s>0Id8fU7;b{x5xNPZ1S6uyUa&`Ru zQ*NHSaN^*>gI$Yr%KadN2YyL>0^j_|fopA-9)!nD$hdOxZWA+K$chAF%?)**Hs_A;hUUvqq7EgS4p6PiG&^0#*!%r2VDQ5F1b(_|+p^)}R87}0_QGj)x{k4@ z{IYv^&Caar=%ZVAw|=QF;?l=9ZJBrTGhNH3a?I}N*4+oL2mgP~WdC#f&5mzO%q?DH zGPf#cZdL9&#$MKS^zr-H-P!lD%Cq0s(U0%p>~z1Uu{W8Knr>^G@y(lm^_rUku z@5(Q`qgQtir`hQ`#-8%a?%^~$v!0_LUTd=GJWF_Wr|TGd=+@ntUzNL_v6poneRS*Y z%&*F`uA?7bV>(7^=*&9n)-|K$dsR)>G4{}_ySt9Dr~I-zdUf}3nw_p=>?yzO9!}Y; z!J_jFH(k8(!acw2S~itqc1O4F?s~=^x^;KvSLLo}>}6d?AKj{*xn;RyrH0OH*N)A) zs(6jb+^U?pRk`aJd+62OUB}o{e%T$px_dayPUV>0Yq01%!+Y^+?4etCXMS1kI>uhsbM(=zW2J^pcy*`i7<=Pk!)weM3_8zn=I-k^YEjhDzWz%(xJ@o4Cu4C*mw<>3DRjwSfJGymu*E9B*U(9Q;=sd$)=C7G^_fuVq zllfISx^;KgGxpG}+Ld2+SAN-@xmCIA7<=f|-5o16bpG@5Rr?n`k5}}|cIB7dUB}o% zukNlKvpe&va^_d%$}zj6TXzr7HCS|>A-rZM^Q&^^m*wcy-Cf7nV{TQh{Ia|9%kIpr z%F(f7q^8bAH?A1p`-`qsGwV6}$}hV!zbbbC#%&*FsUzIb*HCdeOS^oE~ zD^{*7PSsS7@g7dIlety7^2_e%)!kjk*kf*0uKcpQ@*DTbjMQ|yWZ%;Brx&j{m0xy8 zx9*;G9ew4P-I-sNyPmO!ZrxouW_RY7{-75iY3FhJbAQN^vibUR^_f^?1j_p zWNuZi{IWZGb$8b>_Ly6hhu_JJ)N~84*~$E>TsdZUbnEU}*U?vw*`4`Sx$8N4=+)hs z-?4bA z*gDa*YAVO<9$vGP`Bk}c%~#Uv_71Ri5=6edbr?%&*E_&)Cbljy^h`j7Uwl zbrWyhv#NZ*x>Gr3cXX?E=2qp(FT01+>{Nc)UHN5q^y=>66lYBqoo6|5_Boqo6|Zc% zjT}Pk!RXOvka^=}EQbXq-N9TR`3r}<{o6IlE(W|>FzwFN3 zs+_r1x$?{I;WRsyUv^i1@m_;P=NVqNdhXe07q2myTa`1vDtA3&kNH)(a?I|^F}pLr zD$lx(KJ%+`$4U*IaGIUUFS{$h?2caDUHN5q=2qp*t;&^Ob`Pi7iGFLa=sd%wO{YJ- z?eVU4XV!J}nO~JNzbbbtcyNbyvAg1Rj&N9yXzQx;WRsyUv^i1*&V&QJM*h@ z=2zu4SahCY+r`sPI_U_n;ntnXF}pLrDo3~Oo^>64=9lHnt;$`;*i(Ml9sN31YUqSl zce;+T$K0x1`DJ(4G4{e~b~3*zSB}{o-MTyTTZ2XC8M?IxAK3V_u2s|Zj6HPg?#eOV znOl{kS9i~Pj=u8C?#!*qm1D<94V`a){i(5Yic>a~Uv@{Y?(RCq9=dgR<(S=F&)5sE z*~$E>Tse+=4HlhexODepmu>%P*WzS;SZcX&e%Wq(&zDte%a3as$4l{cXaFSu4n9_TX$EE+1>Swz3`fy;kO2h&NGBpcQUss zSAN-@xmCIA7<=f|-Cf7nV{TQBUfrE}cC6IUxnlm_0|yRtt(vZ9>?z0W&itxeIc9fs z>+Y^+?4etCSB}|Td9J~t^9)a1zUz)P#Vbzbm)+5;yEC^cSAN-@xmCIA7<=f|-Cf7n zV{S2bjMUIsed7bGPAN`tqFZ-oepT*z#-4J_?#!>sm1A~Cx9;wGjvjh-_ZloZ&oI2@ zf7(qyIod1wWxMMbd&)1nqgQulZdI=QvO9CDa@R5T!fSShXU9qno$#8S%&*GPt-CY7 zDtA3&PdUarbE|UYm)+5;yED%z5R#IGJ0OGq)5GErq6CT%}&=b_Ly6hGq)3D zRjxc|%=FoP&+6-*xvY3)({+r!aGIU2W9%`vDrat0p7k7kbnEWQF}o|raqrV$=Di=> z^x3Twk9I9i=2zwD*4?wNqmORYuKcpQ^2_ewG&^0#*kf*4p0Uzrcgg*qShub?jfq~} zUHN5q*D?0OX?7~d?9TkEocUFG)^+sJt-JSWF!SC&Y*voG`UtP#R!!G2_Ly6hqgQv& zdX7GGt8(R+-IZT<52x9Qjx%QZ?1oo&D#z@OZrxouW_Q;!_QGp+D#z^3{HmP!RXKCr zr@_p7|J_;3Pr2=pu4U78j6LR7<*sAwF}Es5ukN1p9DU|i<;pL+E5C8inCY|o=2zwD*4>q3c1O4Ft{k(w>lu6DH9M7Kc4vNB-lsw5y+5(`vM=oTVb{h)ukNn= zvb*aTd(5rMUB}pCZdHz6-976$`pm7$9V<0-UUvKYKls5TyoOtMx}LGe{Hh$?y1R1B z?&#Ltm1A~yJ!3DtW@k9A!J_jF;nkh4W9*?p0UULsvOC#=vM8@FT1;rYq01%!znXIAHS}6 z#i{(VJ9>3@*D?0ctGg?|?Cv_o9&@X5<(S=>UzR&o`s|F&zwpjei&LEF*4>$3mAjs? zhi=vGI>sJzt8(<}?#eH_Gq)=5(_rSk|7>ad$??5Mdqux&cO7F-`DJ(X>h8)hyQ5oo zSB}};^^85{SLNY2W2Vn;cy%Xpt8(<}?#!*qUB}pCZdL9&#vXI4a`fu%%x|9tGw=Pc zH=cXJ1?BryP1iH_lw)>xJ!4NfW_NV!?#eN{qg!`Zj@ey#&Y0=5d-j&UnEv4ZcP&ok zm)+5;yEC^cN3ZVA+^XDlj6LR7<*sAwF}E!5(_rSkH}8Aneanj1nCRBsnO~K=p0TGK zv%BjVd&)7pqg!`Zj@cdEx_ib?z0Wj&9u@9lsDeef*jK zT6pTEH-dC$LIJ~DfD@ye!h z%w&z%kJpa-I-gJhu7@v!(Qe+Kd}9eZ@aKKjTv6EQ#odL<(S>Wt2WYj$Q`N1ypsIrFP> z*E9Cet-FWkJ`FnW9Zs{;b&NgcR^`mC%CnxMA5ODV`DJ(Im)+5;yQ5>rN)4S8zVgU* zbISLsJF~8%k8a&vIc9g|nBBu`c4l2ipZQfe^Q&^^xdw~QGu&|h!;9|xzpiD|b&S1m znw_p=>@l}0XKq!V^&I_hnw`ooyDPtOcdXRV85!C+a#8V$llfISx^?%g>*%9fcUO+t zT{&j=@S2@j*U@KwSzd!h=NUeK){b4fb{*{%{jy#8Wp~#x_QGj)x{k5O+^U?pRe9EP z^uuX(hF`}@pPfIybRnA;z%=Fp4ZSA&Y!{z&xP1iB@!fAFg zw<>3DRnFY1JnK37=+)hoUv^i1;}!j~UHN5q*D?0OX?8NVDrat0&fKaz>pA-9 z)!oByp9Y=x{>z!;n>T;2Yr|`HGQTQkepQZc-976%`pPl8E641v9J70P%}(^}nCY_< zUfr4X9DU`N-IZT3YUqc+F1c zSLMvF%F(U6XI)2MIc9g|nBA4*jF~>W?_Ij_!aX~?7N_#d?&#Isv!0``{Ia|9%kHjY z?1j_pWNuZ?+_Jn+gPHe!=*F>GR~4@@(XG2Pzbbbuf&%}#XNr$OhvFW++0j%nq4H9NDe zqmORY&fKcpb&S1mnw`w8%9&f0qgQulo*gSSboTGNe#3_EcCDIO&(RO3*{K|}yK>C# z=+@n{uA{FUv%7N4?#gow7M*8UF|uySnZ+wk<(J*VX?A8kN1wS>IdiLW*D?0OX?8NV zDrat4?pUdz^Q(nx=G^_Au8oOq-JSVWdDeCG!)tab$Ly{gvpc$V_pIyaE640!gGJ{V z!fAH8j9h0c<*W8DD&MQwnROj~bnEWS zugbHoqaR+gQ#odL<(S>kt-CYVeOh$hd-^Ra#`iwhwQ6QPM?aiqr|TGd$}hWz)9lQ8 zjy`j%a^_a$%CTdmhR*u!SFBuFyt1hrvwL{W&aCU`qg!`pepQ}z9sTf{oysx0E5~uK z!J_jFuY7Xp`O}Nnn9QxpnOl`-Jx4#BW~b{Id&)1nhturLdX7GGt8&Lm4V^E}UNU^k z1J8L4x9(Jq*@l}0XKq!V z^&I_hnw_p=>?yzO9!|59d7jKl4YxDivv~Ww^1Z5N)^qgHtGg?|?5_N>dpOO`tmo*X zS9fP_RjxeOV9|Moe_OrihOOW3S~itqb`P)F>3YT<^Q&^^SLIpP(GRcL>3YVVavXQZ zN{t=nR^`mC%CnxMk6ztf`DJ(Im)*l@c4j?CAHBMJjTVo^GbzXHt{k&_c+F1NGxnHY zl{3F8&$^C&c+Jl6d@O8C<-e1;RXKC3@~r3RqgQuVe%W35W%qEJoy_x8wpGKkE641v z9J70P%}&=d_LyIlGruZVu50+6zZ9N)^cypMd)bOm(ln1>m=<%>V!Z literal 82345 zcmeI!F^e2^7{~Ev?~d~zXNW-!LI@jD84HUvK}5Vng@A}C2!b|33P~E1_KFZAXQP6p zr1J&HktzhiS0Gr0KwgD~jR^kt{0B107@i+-_0yQs(CmgBW@{9L*Hz5K4%$5+qWeC+sz#YP+{kNsJGU#gFfcu5=f_h||C~Jfpgi}=GUA`NB7XQh;`DCBou4A!y&G}=e#Do5 zMg0BZBECCUPkM0X(xmb8&57Ih-v9L4+mCMRlm7SS_V#PnFTZze=jT5@`|85(8#|Al zxH@aBu6d)rL%m(w*xx%|{<7HG-`kAiv00Y1i0{wD_UiW6BO3>|ub%l}=biKCmai`_ zTwLy4{o?Yq>$?Y!lxuST_vU{CHgDXAX1);=Xp)lejo=CyahL+KCHV!zTo>;;SOHl` zk1rT30&mg=Ysu?bz#sPRo`4|Bd0LxXc(#}uyJkW?LM3U%jv4TiGj zyp9Q7JF;Rz+jP1VWFaehNr0?gvO0Zv3^WU}kafu4S&)L%A(N_Kj>tk5vdX7mp#LoJ zwt9KKv-D;-%QQn)rV^HfN`|EBNUTvU35@I1M6U&fP9f|#A_UBv4 zP?m2jCKIqDONfOmAr`WPSjZA$AxnsbEFl)MgjmSR#DY|&5|;Ez$kLnPEYl2GnMzpF zDU35@I1M z6AMzAN?6h>Axm$DvrIE&Wh!AwuY@eU8O}1zkd>)~C7}|egjmQDVj)Y2g)AW!vV>U3 z5@I1sh=nX67P5p`$P!{X$zpNB;(TK{k*c%!#$qzTJhFsX$P!{9ONfOmAr`WPSjZA$ zAxnsbtV}FOWh!AwuY@eU8O}1zkd>)~CA|`|^kz8AG(%RV5|;EzPO{jazCXiePNeEA zhLspe$Brz$8O}1zkd>)~CA|`|^kz8AG;@;0#xfhTujE9k&SKw;0W$(knU1Vt@Mne48O_O53z>=0vK_V$p>mx_;D<)wLi7+Y2JA zy$fEtc0pFxf*5Qs=p>6BweM(TU7SeOSw>cnaoax#$ZGF`m#$rq)wLi7+Y2JAy$fEt zc5#x$HoI=qsUWhZv`q_HPNeEA7F`&k>qiY)T?=Bcy&$sMyWpj37bjV4v+FjU3LsFY&DU{{0-Xy~e4WqJ_X++5JfxUd diff --git a/api/app/tests/geospatial/4326_lats.tif b/api/app/tests/geospatial/4326_lats.tif index d11d5347dd2fc0974bbf2f9138c0acbbd0700e6f..4a629497fa1d385e48b27fade642d3709708807a 100644 GIT binary patch literal 80426 zcmeI4|4-C)9LGP$86iT0hThUdy`qsqtZWP`>)o)CvG+Orve4X^4s4VrC5r{=!W396 z8En#M$Fy)#tW1<88;M0@lSbvjFG(U|Uz1GIkE$ z^T;{KALM!exc4j>JG+eotbi8Ea$zY5l)DltK{5R z8Z!FpKbCXwVP@#R2Dvvd({W}S$9W~oaT+$sb*)?-v2!oH^*)rT(#=!d8vVskr8X|yw^2q@ovYNARgtHiW;H(2+{0v&PBPx8cXG_ z6FPdm?SYA5drO1&Cyml^%GTvqPv7-M@bk>>rJ2vZn6_kQ%HwZ4cfJSS_KgQT1iQPl zgKyRkwO7eC!~;9f9X#(pUmhoYKs>Mm-NBK=XEK^^Sp)IF4s-{Zfs>we=>y_{9q0~v zXLa@T{A&%w13SfUGjV&!Vt$}!82fBmE#D-#@^a1g}4s-|Gm)E8B{bLQp13S()R#umjz}x7me*8PW&D13SvmqX z2I7Go=niHS#m}89eLy_01KmM={j3}9SFC|}UMm-NDm|k?EVH4~PeLpgU+^9vU4jeLy_01Koiw8983}r!^1{>_B&LY~OW%vh)G* zzz%c=D|;`m?zm(P!~;9f9rRBcj9)2zKs>Mm-ND+F{-L2EYakxjf$rdHc285W^a1g} z4s-`amB$NXqz{M(cAz_OxBh7A|F8z)fgR`$-Wxa^`h@fW@xTsr2lsl~Iy-;22I7Go z=ng(iZK=wVJ|G_0f$m^dUQ_0kL2Do$*n#fgQ~&3&Tcr<(2X>%4c(h}8zeoCjcwh&* zgXX}FhK2!aARgF(?!b6!^5#n)5D)A?cW~HOm3-u)H4qQ%KzFdHsC;;}^a1g}4s-`6 z>bJB7E?5Kczz%c=Sw}aPl}R5E59~m9aLy_4O_M$#9@v5I;I+iOm_6sMfp}mCx`S)W zbGwc70r9{NbO#$tR`2XOXAQ&yJJ226wg06xYorf|2X>%4c)xdP!j0dofp}mCx`UVq z-<5Ld1LA=l=nniTY0a_H2gCzA&>cLOV>UJRTLbaH4s-{fRnA+mK>C1qUMm-N8dEI_vVJ4~PeLpgY)8(w^&GeLy_01KmM(>-Mft=>y_{9q0}Q2e#GJoU{hwfgR`$ zHh9Xi5~L4^2X>%4xS3iUcj#wpARgF(?qG9X;b4aJ0r9{NbO+(J8yfphSOf9E4s-`q z9cv2;qz{M(cAz`B?^@3MaOnf$fgR`$KJsR{YmZw4@xTsr2M;gG>`0P6ARgF(?qF|G PMn&5(Yakxjf$rcx3V`Vi literal 70478 zcmeIyzf06{90u^ushFjNLKZ}sP*XuPG!%qc4dPIPhBdSVO+`yIMFjo4riO;-4~W(j zA{rc;`U|3=!7+k1p)IiYb6#LUQ|rs)xKCfc_lcJxkbXJ59q^scq*pIW=O-rDo={^K?uT0&c2bE|7_ z?_aoa`pD$wXaA;uek{D7di?$B)ek4n{n*&pI6wRQ#yqp-mk*VF{i5RR{Xh9Zo%VIS zGqTc+%W&`C&kk-MS?=aLqrK18p5DCuV*jgKdEV^niPcar|yORTt{d2LdTj z1X~aRTYx|c6d?k(0D%-JLIi990x3|02-pGyQlJPCumuRDKoKHf3lK~5wHaaq(BiOU<(jPfg(h} z79fxUMTmecKp+K*5CL0&KnfHg0=58w6evOjYyko(P=pBB0t8Z^2obOa2&6y}B47&; zNP!|mz!o5o0!4^`EkGaziVy)?fItcqAp*7lffOi01Z)8UDNuw6*a8Gnpa>DL1qh@- z5h7p<5J-U{M8FmxkOD=BfGt2E1&R;>TYx|c6d?k(0D%-JLIi990x3|02-pGyQlJPC zumuRDKoKHf3lK~5wHaaq(BiOU<(jPfg(h}79fxUMTmecKp+K*5CL0&KnfHg0=58w6evOjYyko( zP=pBB0t8Z^2obOa2&6y}B47&;NP!|mz!o5o0!4^`EkGaziVy)?fItcqAp*7lffOi0 z1Z)8UDNuw6*a8Gnpa>DL1qh@-5h7p<5J-U{M8FmxkOD=BfGt2E1&R;>TYx|c6d?k( P0D%-JLIiBVe-r!yJrQ4s diff --git a/api/app/tests/geospatial/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py index 55fc8f134..811101517 100644 --- a/api/app/tests/geospatial/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -173,6 +173,8 @@ def test_latitude_array(): lats_4326_tif = os.path.join(os.path.dirname(__file__), "4326_lats.tif") with WPSDataset(ds_path=lats_3005_tif) as lats_3005_ds, WPSDataset(ds_path=lats_4326_tif) as lats_4326_ds: output_ds: WPSDataset = lats_3005_ds.warp_to_match(lats_4326_ds, "/vsimem/test_lats.tif") - original_lats = lats_4326_ds.generate_latitude_array() + original_ds = gdal.Open(lats_4326_tif) + original_lats = original_ds.GetRasterBand(1).ReadAsArray() warped_lats = output_ds.generate_latitude_array() assert np.all(original_lats == warped_lats) == True + output_ds = None From 917845153b5523d2d470b823c84d056cb7c9a782 Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Mon, 21 Oct 2024 14:56:33 -0700 Subject: [PATCH 15/16] Rename arguments for clarity --- api/app/geospatial/wps_dataset.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/api/app/geospatial/wps_dataset.py b/api/app/geospatial/wps_dataset.py index 62e931704..bff7ebcf0 100644 --- a/api/app/geospatial/wps_dataset.py +++ b/api/app/geospatial/wps_dataset.py @@ -97,18 +97,18 @@ def warp_to_match(self, other, output_path: str, resample_method: GDALResampling """ Warp the dataset to match the extent, pixel size, and projection of the other dataset. - :param ds_to_match: the reference dataset raster to match the source against + :param other: the reference WPSDataset raster to match the source against :param output_path: output path of the resulting raster :param resample_method: gdal resampling algorithm :return: warped raster dataset """ - source_geotransform = other.ds.GetGeoTransform() - x_res = source_geotransform[1] - y_res = -source_geotransform[5] - minx = source_geotransform[0] - maxy = source_geotransform[3] - maxx = minx + source_geotransform[1] * other.ds.RasterXSize - miny = maxy + source_geotransform[5] * other.ds.RasterYSize + dest_geotransform = other.ds.GetGeoTransform() + x_res = dest_geotransform[1] + y_res = -dest_geotransform[5] + minx = dest_geotransform[0] + maxy = dest_geotransform[3] + maxx = minx + dest_geotransform[1] * other.ds.RasterXSize + miny = maxy + dest_geotransform[5] * other.ds.RasterYSize extent = [minx, miny, maxx, maxy] # Warp to match input option parameters From 2cd62a92dde4eb0b17086206d48afa927696396f Mon Sep 17 00:00:00 2001 From: Conor Brady Date: Mon, 21 Oct 2024 16:00:47 -0700 Subject: [PATCH 16/16] update replace no data value impl and test --- api/app/geospatial/wps_dataset.py | 16 +++-- api/app/tests/geospatial/test_wps_dataset.py | 62 +++++++++++++++----- 2 files changed, 54 insertions(+), 24 deletions(-) diff --git a/api/app/geospatial/wps_dataset.py b/api/app/geospatial/wps_dataset.py index bff7ebcf0..1e46e0e42 100644 --- a/api/app/geospatial/wps_dataset.py +++ b/api/app/geospatial/wps_dataset.py @@ -125,21 +125,19 @@ def warp_to_match(self, other, output_path: str, resample_method: GDALResampling ) return WPSDataset(ds_path=None, ds=warped_ds) - def replace_nodata_with(self, new_no_data_value: int): + def replace_nodata_with(self, new_no_data_value: int = 0): """ - Modifies the dataset inplace by replacing the nodata value with the supplied one - + Reads the first band of a dataset, replaces NoData values with new_no_data_value, returns the array and the nodata value. :param new_no_data_value: the new nodata value """ - band: gdal.Band = self.ds.GetRasterBand(self.band) + + band: gdal.Band = self.ds.GetRasterBand(1) nodata_value = band.GetNoDataValue() array = band.ReadAsArray() - if nodata_value is not None: - modified_array = np.where(array != nodata_value, array, new_no_data_value) - band.WriteArray(modified_array) - band.SetNoDataValue(new_no_data_value) - band.FlushCache() + array[array == nodata_value] = new_no_data_value + + return array, new_no_data_value def generate_latitude_array(self): """ diff --git a/api/app/tests/geospatial/test_wps_dataset.py b/api/app/tests/geospatial/test_wps_dataset.py index 811101517..75d20f80a 100644 --- a/api/app/tests/geospatial/test_wps_dataset.py +++ b/api/app/tests/geospatial/test_wps_dataset.py @@ -10,7 +10,7 @@ zero_tif = os.path.join(os.path.dirname(__file__), "zero_layer.tif") -def create_test_dataset(filename, width, height, extent, projection, data_type=gdal.GDT_Float32, fill_value=None) -> gdal.Dataset: +def create_test_dataset(filename, width, height, extent, projection, data_type=gdal.GDT_Float32, fill_value=None, no_data_value=None) -> gdal.Dataset: """ Create a test GDAL dataset. """ @@ -32,20 +32,44 @@ def create_test_dataset(filename, width, height, extent, projection, data_type=g # Create some test data (e.g., random values) rng = np.random.default_rng(seed=42) # Reproducible random generator - random_data = rng.random((height, width)).astype(np.float32) + fill_data = rng.random((height, width)).astype(np.float32) - if fill_value is None: - # Write data to the dataset - dataset.GetRasterBand(1).WriteArray(random_data) - else: - fill_data = np.full_like(random_data, fill_value) - dataset.GetRasterBand(1).WriteArray(fill_data) + if fill_value is not None: + fill_data = np.full((height, width), fill_value) dataset.GetRasterBand(1).SetNoDataValue(0) + dataset.GetRasterBand(1).WriteArray(fill_data) return dataset +# def create_test_dataset_with_no_data_value(filename, width, height, extent, projection, data_type=gdal.GDT_Float32, fill_value: int) -> gdal.Dataset: +# """ +# Create a test GDAL dataset. +# """ +# # Create a new GDAL dataset +# driver: gdal.Driver = gdal.GetDriverByName("MEM") +# dataset: gdal.Dataset = driver.Create(filename, width, height, 1, data_type) + +# # Set the geotransform +# xmin, xmax, ymin, ymax = extent +# xres = (xmax - xmin) / width +# yres = (ymax - ymin) / height +# geotransform = (xmin, xres, 0, ymax, 0, -yres) # Top-left corner +# dataset.SetGeoTransform(geotransform) + +# # Set the projection +# srs = osr.SpatialReference() +# srs.ImportFromEPSG(projection) +# dataset.SetProjection(srs.ExportToWkt()) + +# rng = np.random.default_rng(seed=42) # Reproducible random generator +# random_data = rng.random((height, width)).astype(np.float32) + +# fill_data = np.full_like(random_data, fill_value) +# dataset.GetRasterBand(1).WriteArray(fill_data) + + def test_raster_with_context(): """ with opens the dataset and closes after the context ends @@ -57,13 +81,21 @@ def test_raster_with_context(): def test_raster_set_no_data_value(): - extent = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax - ds_1 = create_test_dataset("test_dataset_no_data_value.tif", 1, 1, extent, 4326, data_type=gdal.GDT_Byte, fill_value=2) - with WPSDataset(ds_path=None, ds=ds_1) as wps_ds: - assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == 0 - - wps_ds.replace_nodata_with(-1) - assert wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() == -1 + original_no_data_value = 0 + driver: gdal.Driver = gdal.GetDriverByName("MEM") + dataset: gdal.Dataset = driver.Create("test_dataset_no_data_value.tif", 2, 2, 1, eType=gdal.GDT_Int32) + fill_data = np.full((2, 2), 2) + fill_data[0, 0] = original_no_data_value + dataset.GetRasterBand(1).SetNoDataValue(original_no_data_value) + dataset.GetRasterBand(1).WriteArray(fill_data) + + with WPSDataset(ds_path=None, ds=dataset) as wps_ds: + original_array = wps_ds.as_gdal_ds().GetRasterBand(1).ReadAsArray() + original_nodata_value = wps_ds.as_gdal_ds().GetRasterBand(1).GetNoDataValue() + updated_array, updated_nodata_value = wps_ds.replace_nodata_with(-1) + + assert original_array[0, 0] == original_nodata_value + assert updated_array[0, 0] == updated_nodata_value def test_raster_mul():