From 6c664307ba7ca34422da082f71e09a74b677a1fe Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Wed, 28 Aug 2024 09:51:57 -0700 Subject: [PATCH 01/32] Moved scale_meters argument into class definition --- city_metrix/layers/albedo.py | 10 +- city_metrix/layers/alos_dsm.py | 5 +- .../layers/average_net_building_height.py | 8 +- city_metrix/layers/built_up_height.py | 5 +- city_metrix/layers/esa_world_cover.py | 24 ++-- .../layers/land_surface_temperature.py | 6 +- city_metrix/layers/landsat_collection_2.py | 5 +- city_metrix/layers/nasa_dem.py | 5 +- city_metrix/layers/natural_areas.py | 6 +- city_metrix/layers/ndvi_sentinel2_gee.py | 5 +- city_metrix/layers/sentinel_2_level_2.py | 2 - city_metrix/layers/tree_canopy_height.py | 7 +- city_metrix/layers/tree_cover.py | 6 +- city_metrix/layers/urban_land_use.py | 5 +- city_metrix/layers/world_pop.py | 6 +- tests/test_layer_dimensions.py | 1 + tests/test_layer_parameters.py | 103 ++++++++++++++++++ 17 files changed, 161 insertions(+), 48 deletions(-) create mode 100644 tests/test_layer_parameters.py diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py index 7bf7b11..9297230 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -4,12 +4,12 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection - class Albedo(Layer): - def __init__(self, start_date="2021-01-01", end_date="2022-01-01", threshold=None, **kwargs): + def __init__(self, start_date="2021-01-01", end_date="2022-01-01", scale_meters=10, threshold=None, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date + self.scale_meters = scale_meters self.threshold = threshold def get_data(self, bbox): @@ -30,7 +30,7 @@ def mask_and_count_clouds(s2wc, geom): nb_cloudy_pixels = is_cloud.reduceRegion( reducer=ee.Reducer.sum().unweighted(), geometry=geom, - scale=10, + scale=self.scale_meters, maxPixels=1e9 ) return s2wc.updateMask(is_cloud.eq(0)).set('nb_cloudy_pixels', @@ -88,7 +88,9 @@ def calc_s2_albedo(image): s2_albedo = dataset.map(calc_s2_albedo) albedo_mean = s2_albedo.reduce(ee.Reducer.mean()) - data = get_image_collection(ee.ImageCollection(albedo_mean), bbox, 10, "albedo").albedo_mean + data = (get_image_collection( + ee.ImageCollection(albedo_mean), bbox, self.scale_meters, "albedo") + .albedo_mean) if self.threshold is not None: return data.where(data < self.threshold) diff --git a/city_metrix/layers/alos_dsm.py b/city_metrix/layers/alos_dsm.py index 95187e6..024a9f8 100644 --- a/city_metrix/layers/alos_dsm.py +++ b/city_metrix/layers/alos_dsm.py @@ -6,8 +6,9 @@ class AlosDSM(Layer): - def __init__(self, **kwargs): + def __init__(self, scale_meters=30, **kwargs): super().__init__(**kwargs) + self.scale_meters = scale_meters def get_data(self, bbox): dataset = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2") @@ -16,6 +17,6 @@ def get_data(self, bbox): .select('DSM') .mean() ) - data = get_image_collection(alos_dsm, bbox, 30, "ALOS DSM").DSM + data = get_image_collection(alos_dsm, bbox, self.scale_meters, "ALOS DSM").DSM return data diff --git a/city_metrix/layers/average_net_building_height.py b/city_metrix/layers/average_net_building_height.py index a5b26f0..2a215c7 100644 --- a/city_metrix/layers/average_net_building_height.py +++ b/city_metrix/layers/average_net_building_height.py @@ -5,10 +5,10 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection - class AverageNetBuildingHeight(Layer): - def __init__(self, **kwargs): + def __init__(self, scale_meters=100, **kwargs): super().__init__(**kwargs) + self.scale_meters = scale_meters def get_data(self, bbox): # https://ghsl.jrc.ec.europa.eu/ghs_buH2023.php @@ -18,6 +18,8 @@ def get_data(self, bbox): # GLOBE - ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A") anbh = ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A") - data = get_image_collection(ee.ImageCollection(anbh), bbox, 100, "average net building height").b1 + data = (get_image_collection( + ee.ImageCollection(anbh), bbox, self.scale_meters, "average net building height") + .b1) return data diff --git a/city_metrix/layers/built_up_height.py b/city_metrix/layers/built_up_height.py index 05f3159..bc89d70 100644 --- a/city_metrix/layers/built_up_height.py +++ b/city_metrix/layers/built_up_height.py @@ -7,8 +7,9 @@ class BuiltUpHeight(Layer): - def __init__(self, **kwargs): + def __init__(self, scale_meters=100, **kwargs): super().__init__(**kwargs) + self.scale_meters = scale_meters def get_data(self, bbox): # Notes for Heat project: @@ -18,6 +19,6 @@ def get_data(self, bbox): # ee.ImageCollection("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_R2023A") built_height = ee.Image("JRC/GHSL/P2023A/GHS_BUILT_H/2018") - data = get_image_collection(ee.ImageCollection(built_height), bbox, 100, "built up height") + data = get_image_collection(ee.ImageCollection(built_height), bbox, self.scale_meters, "built up height") return data.built_height diff --git a/city_metrix/layers/esa_world_cover.py b/city_metrix/layers/esa_world_cover.py index ed5f823..b306d8c 100644 --- a/city_metrix/layers/esa_world_cover.py +++ b/city_metrix/layers/esa_world_cover.py @@ -19,32 +19,32 @@ class EsaWorldCoverClass(Enum): MANGROVES = 95 MOSS_AND_LICHEN = 100 - class EsaWorldCover(Layer): STAC_CATALOG_URI = "https://services.terrascope.be/stac/" STAC_COLLECTION_ID = "urn:eop:VITO:ESA_WorldCover_10m_2020_AWS_V1" STAC_ASSET_ID = "ESA_WORLDCOVER_10M_MAP" - def __init__(self, land_cover_class=None, year=2020, **kwargs): + def __init__(self, land_cover_class=None, year=2020, scale_meters=10, **kwargs): super().__init__(**kwargs) self.land_cover_class = land_cover_class self.year = year + self.scale_meters = scale_meters def get_data(self, bbox): if self.year == 2020: data = get_image_collection( - ee.ImageCollection("ESA/WorldCover/v100"), - bbox, - 10, - "ESA world cover" - ).Map + ee.ImageCollection("ESA/WorldCover/v100"), + bbox, + self.scale_meters, + "ESA world cover" + ).Map elif self.year == 2021: data = get_image_collection( - ee.ImageCollection("ESA/WorldCover/v200"), - bbox, - 10, - "ESA world cover" - ).Map + ee.ImageCollection("ESA/WorldCover/v200"), + bbox, + self.scale_meters, + "ESA world cover" + ).Map if self.land_cover_class: data = data.where(data == self.land_cover_class.value) diff --git a/city_metrix/layers/land_surface_temperature.py b/city_metrix/layers/land_surface_temperature.py index 7fff632..da4922c 100644 --- a/city_metrix/layers/land_surface_temperature.py +++ b/city_metrix/layers/land_surface_temperature.py @@ -5,12 +5,12 @@ import ee import xarray - class LandSurfaceTemperature(Layer): - def __init__(self, start_date="2013-01-01", end_date="2023-01-01", **kwargs): + def __init__(self, start_date="2013-01-01", end_date="2023-01-01", scale_meters=30, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date + self.scale_meters = scale_meters def get_data(self, bbox): def cloud_mask(image): @@ -31,5 +31,5 @@ def apply_scale_factors(image): .map(apply_scale_factors) \ .reduce(ee.Reducer.mean()) - data = get_image_collection(ee.ImageCollection(l8_st), bbox, 30, "LST").ST_B10_mean + data = get_image_collection(ee.ImageCollection(l8_st), bbox, self.scale_meters, "LST").ST_B10_mean return data diff --git a/city_metrix/layers/landsat_collection_2.py b/city_metrix/layers/landsat_collection_2.py index 248227a..d82180d 100644 --- a/city_metrix/layers/landsat_collection_2.py +++ b/city_metrix/layers/landsat_collection_2.py @@ -29,7 +29,8 @@ def get_data(self, bbox): fail_on_error=False, ) - # TODO: Determine how to output xarray - qa_lst = lc2.where((lc2.qa_pixel & 24) == 0) return qa_lst.drop_vars("qa_pixel") + + + diff --git a/city_metrix/layers/nasa_dem.py b/city_metrix/layers/nasa_dem.py index 49b2d13..35daa0c 100644 --- a/city_metrix/layers/nasa_dem.py +++ b/city_metrix/layers/nasa_dem.py @@ -6,8 +6,9 @@ class NasaDEM(Layer): - def __init__(self, **kwargs): + def __init__(self, scale_meters=30, **kwargs): super().__init__(**kwargs) + self.scale_meters = scale_meters def get_data(self, bbox): dataset = ee.Image("NASA/NASADEM_HGT/001") @@ -16,6 +17,6 @@ def get_data(self, bbox): .select('elevation') .mean() ) - data = get_image_collection(nasa_dem, bbox, 30, "NASA DEM").elevation + data = get_image_collection(nasa_dem, bbox, self.scale_meters, "NASA DEM").elevation return data diff --git a/city_metrix/layers/natural_areas.py b/city_metrix/layers/natural_areas.py index 1dfbc73..61c648a 100644 --- a/city_metrix/layers/natural_areas.py +++ b/city_metrix/layers/natural_areas.py @@ -4,13 +4,13 @@ from .layer import Layer from .esa_world_cover import EsaWorldCover, EsaWorldCoverClass - class NaturalAreas(Layer): - def __init__(self, **kwargs): + def __init__(self, scale_meters=10, **kwargs): super().__init__(**kwargs) + self.scale_meters = scale_meters def get_data(self, bbox): - esa_world_cover = EsaWorldCover().get_data(bbox) + esa_world_cover = EsaWorldCover(scale_meters=self.scale_meters).get_data(bbox) reclass_map = { EsaWorldCoverClass.TREE_COVER.value: 1, EsaWorldCoverClass.SHRUBLAND.value: 1, diff --git a/city_metrix/layers/ndvi_sentinel2_gee.py b/city_metrix/layers/ndvi_sentinel2_gee.py index c5b21b9..ca7fc05 100644 --- a/city_metrix/layers/ndvi_sentinel2_gee.py +++ b/city_metrix/layers/ndvi_sentinel2_gee.py @@ -10,9 +10,10 @@ class NdviSentinel2(Layer): Notebook: https://github.com/wri/cities-cities4forests-indicators/blob/dev-eric/scripts/extract-VegetationCover.ipynb Reference: https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index """ - def __init__(self, year=None, **kwargs): + def __init__(self, year=None, scale_meters=10, **kwargs): super().__init__(**kwargs) self.year = year + self.scale_meters = scale_meters def get_data(self, bbox): if self.year is None: @@ -39,7 +40,7 @@ def calculate_ndvi(image): ndvi_mosaic = ndvi.qualityMosaic('NDVI') ic = ee.ImageCollection(ndvi_mosaic) - ndvi_data = get_image_collection(ic, bbox, 10, "NDVI") + ndvi_data = get_image_collection(ic, bbox, self.scale_meters, "NDVI") xdata = ndvi_data.to_dataarray() diff --git a/city_metrix/layers/sentinel_2_level_2.py b/city_metrix/layers/sentinel_2_level_2.py index a609293..a7ae944 100644 --- a/city_metrix/layers/sentinel_2_level_2.py +++ b/city_metrix/layers/sentinel_2_level_2.py @@ -50,6 +50,4 @@ def get_data(self, bbox): cloud_masked = s2.where(s2 != 0).where(s2.scl != 3).where(s2.scl != 8).where(s2.scl != 9).where( s2.scl != 10) - # TODO: Determine how to output as an xarray - return cloud_masked.drop_vars("scl") diff --git a/city_metrix/layers/tree_canopy_height.py b/city_metrix/layers/tree_canopy_height.py index 2785272..9141f84 100644 --- a/city_metrix/layers/tree_canopy_height.py +++ b/city_metrix/layers/tree_canopy_height.py @@ -5,19 +5,20 @@ import xee import ee - class TreeCanopyHeight(Layer): name = "tree_canopy_height" NO_DATA_VALUE = 0 - def __init__(self, **kwargs): + def __init__(self, scale_meters=1, **kwargs): super().__init__(**kwargs) + self.scale_meters = scale_meters def get_data(self, bbox): canopy_ht = ee.ImageCollection("projects/meta-forest-monitoring-okw37/assets/CanopyHeight") # aggregate time series into a single image canopy_ht = canopy_ht.reduce(ee.Reducer.mean()).rename("cover_code") - data = get_image_collection(ee.ImageCollection(canopy_ht), bbox, 1, "tree canopy height") + data = get_image_collection(ee.ImageCollection(canopy_ht), bbox, + self.scale_meters, "tree canopy height") return data.cover_code diff --git a/city_metrix/layers/tree_cover.py b/city_metrix/layers/tree_cover.py index fadc9bb..9caa483 100644 --- a/city_metrix/layers/tree_cover.py +++ b/city_metrix/layers/tree_cover.py @@ -5,7 +5,6 @@ import xee import ee - class TreeCover(Layer): """ Merged tropical and nontropical tree cover from WRI @@ -13,10 +12,11 @@ class TreeCover(Layer): NO_DATA_VALUE = 255 - def __init__(self, min_tree_cover=None, max_tree_cover=None, **kwargs): + def __init__(self, min_tree_cover=None, max_tree_cover=None, scale_meters=10, **kwargs): super().__init__(**kwargs) self.min_tree_cover = min_tree_cover self.max_tree_cover = max_tree_cover + self.scale_meters = scale_meters def get_data(self, bbox): tropics = ee.ImageCollection('projects/wri-datalab/TropicalTreeCover') @@ -24,7 +24,7 @@ def get_data(self, bbox): merged_ttc = tropics.merge(nontropics) ttc_image = merged_ttc.reduce(ee.Reducer.mean()).rename('ttc') - data = get_image_collection(ee.ImageCollection(ttc_image), bbox, 10, "tree cover").ttc + data = get_image_collection(ee.ImageCollection(ttc_image), bbox, self.scale_meters, "tree cover").ttc if self.min_tree_cover is not None: data = data.where(data >= self.min_tree_cover) diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py index f34408a..da83219 100644 --- a/city_metrix/layers/urban_land_use.py +++ b/city_metrix/layers/urban_land_use.py @@ -7,9 +7,10 @@ class UrbanLandUse(Layer): - def __init__(self, band='lulc', **kwargs): + def __init__(self, band='lulc', scale_meters=5, **kwargs): super().__init__(**kwargs) self.band = band + self.scale_meters = scale_meters def get_data(self, bbox): dataset = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1") @@ -27,6 +28,6 @@ def get_data(self, bbox): .rename('lulc') ) - data = get_image_collection(ulu, bbox, 5, "urban land use").lulc + data = get_image_collection(ulu, bbox, self.scale_meters, "urban land use").lulc return data diff --git a/city_metrix/layers/world_pop.py b/city_metrix/layers/world_pop.py index f3cf1a2..d77cc2e 100644 --- a/city_metrix/layers/world_pop.py +++ b/city_metrix/layers/world_pop.py @@ -5,10 +5,10 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection - class WorldPop(Layer): - def __init__(self, **kwargs): + def __init__(self, scale_meters=100, **kwargs): super().__init__(**kwargs) + self.scale_meters = scale_meters def get_data(self, bbox): # load population @@ -20,5 +20,5 @@ def get_data(self, bbox): .mean() ) - data = get_image_collection(world_pop, bbox, 100, "world pop") + data = get_image_collection(world_pop, bbox, self.scale_meters, "world pop") return data.population diff --git a/tests/test_layer_dimensions.py b/tests/test_layer_dimensions.py index 15768d6..ee32af2 100644 --- a/tests/test_layer_dimensions.py +++ b/tests/test_layer_dimensions.py @@ -16,3 +16,4 @@ def test_ndvi_dimensions(): assert actual_min == expected_min assert actual_max == expected_max + diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py new file mode 100644 index 0000000..bef67ec --- /dev/null +++ b/tests/test_layer_parameters.py @@ -0,0 +1,103 @@ +from city_metrix.layers import ( + Albedo, + AlosDSM, + AverageNetBuildingHeight, + EsaWorldCover, + EsaWorldCoverClass, + LandSurfaceTemperature, + NasaDEM, + NaturalAreas, + TreeCanopyHeight, + TreeCover, + UrbanLandUse, + WorldPop, NdviSentinel2, BuiltUpHeight +) +from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 + +COUNTRY_CODE_FOR_BBOX = 'BRA' +BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 + +def test_albedo_scale(): + doubled_default_resolution = 2 * Albedo().scale_meters + layer = Albedo(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_alos_dsm_scale(): + doubled_default_resolution = 2 * AlosDSM().scale_meters + layer = AlosDSM(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_average_net_building_height_scale(): + doubled_default_resolution = 2 * AverageNetBuildingHeight().scale_meters + layer = AverageNetBuildingHeight(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_built_up_height_scale(): + doubled_default_resolution = 2 * BuiltUpHeight().scale_meters + layer = BuiltUpHeight(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_esa_world_cover_scale(): + doubled_default_resolution = 2 * EsaWorldCover().scale_meters + layer = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP, scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +# TODO +# def test_high_land_surface_temperature_scale(): + +def test_land_surface_temperature_scale(): + doubled_default_resolution = 2 * LandSurfaceTemperature().scale_meters + layer = LandSurfaceTemperature(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_nasa_dem_scale(): + doubled_default_resolution = 2 * NasaDEM().scale_meters + layer = NasaDEM(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_natural_areas_scale(): + doubled_default_resolution = 2 * NaturalAreas().scale_meters + layer = NaturalAreas(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_ndvi_sentinel2_scale(): + doubled_default_resolution = 2 * NdviSentinel2().scale_meters + layer = NdviSentinel2(year=2023, scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +# TODO +# def test_smart_surface_lulc_scale(): + +def test_tree_canopy_height(): + doubled_default_resolution = 2 * TreeCanopyHeight().scale_meters + layer = TreeCanopyHeight(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_tree_cover(): + doubled_default_resolution = 2 * TreeCover().scale_meters + layer = TreeCover(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_urban_land_use(): + doubled_default_resolution = 2 * UrbanLandUse().scale_meters + layer = UrbanLandUse(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + +def test_world_pop(): + doubled_default_resolution = 2 * WorldPop().scale_meters + layer = WorldPop(scale_meters=doubled_default_resolution) + evaluate_layer(layer, doubled_default_resolution) + + +def evaluate_layer(layer, expected_resolution): + data = layer.get_data(BBOX) + est_actual_resolution = get_resolution_estimate(data) + assert expected_resolution == est_actual_resolution + +def get_resolution_estimate(data): + y_cells = data['y'].size - 1 + y_min = data['y'].values.min() + y_max = data['y'].values.max() + y_diff = y_max - y_min + ry = round(y_diff/y_cells) + return ry From 164897467db2a3598c850dabcd42b420072f83f0 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 29 Aug 2024 15:05:13 -0700 Subject: [PATCH 02/32] Renamed class property in compliance with stac standard. Implemented test_layer_parameters in a manner intended to improve ease of code maintenance. --- city_metrix/layers/albedo.py | 16 +- city_metrix/layers/alos_dsm.py | 11 +- .../layers/average_net_building_height.py | 11 +- city_metrix/layers/built_up_height.py | 11 +- city_metrix/layers/esa_world_cover.py | 15 +- .../layers/land_surface_temperature.py | 13 +- city_metrix/layers/nasa_dem.py | 11 +- city_metrix/layers/natural_areas.py | 11 +- city_metrix/layers/ndvi_sentinel2_gee.py | 10 +- city_metrix/layers/tree_canopy_height.py | 11 +- city_metrix/layers/tree_cover.py | 12 +- city_metrix/layers/urban_land_use.py | 12 +- city_metrix/layers/world_pop.py | 11 +- tests/test_layer_parameters.py | 157 ++++++++++-------- 14 files changed, 197 insertions(+), 115 deletions(-) diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py index 9297230..341786b 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -5,11 +5,19 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection class Albedo(Layer): - def __init__(self, start_date="2021-01-01", end_date="2022-01-01", scale_meters=10, threshold=None, **kwargs): + """ + Attributes: + start_date: starting date for data retrieval + end_date: ending date for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + threshold: threshold value for filtering the retrieval + """ + + def __init__(self, start_date="2021-01-01", end_date="2022-01-01", spatial_resolution=10, threshold=None, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution self.threshold = threshold def get_data(self, bbox): @@ -30,7 +38,7 @@ def mask_and_count_clouds(s2wc, geom): nb_cloudy_pixels = is_cloud.reduceRegion( reducer=ee.Reducer.sum().unweighted(), geometry=geom, - scale=self.scale_meters, + scale=self.spatial_resolution, maxPixels=1e9 ) return s2wc.updateMask(is_cloud.eq(0)).set('nb_cloudy_pixels', @@ -89,7 +97,7 @@ def calc_s2_albedo(image): albedo_mean = s2_albedo.reduce(ee.Reducer.mean()) data = (get_image_collection( - ee.ImageCollection(albedo_mean), bbox, self.scale_meters, "albedo") + ee.ImageCollection(albedo_mean), bbox, self.spatial_resolution, "albedo") .albedo_mean) if self.threshold is not None: diff --git a/city_metrix/layers/alos_dsm.py b/city_metrix/layers/alos_dsm.py index 024a9f8..70000eb 100644 --- a/city_metrix/layers/alos_dsm.py +++ b/city_metrix/layers/alos_dsm.py @@ -6,9 +6,14 @@ class AlosDSM(Layer): - def __init__(self, scale_meters=30, **kwargs): + """ + Attributes: + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, spatial_resolution=30, **kwargs): super().__init__(**kwargs) - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): dataset = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2") @@ -17,6 +22,6 @@ def get_data(self, bbox): .select('DSM') .mean() ) - data = get_image_collection(alos_dsm, bbox, self.scale_meters, "ALOS DSM").DSM + data = get_image_collection(alos_dsm, bbox, self.spatial_resolution, "ALOS DSM").DSM return data diff --git a/city_metrix/layers/average_net_building_height.py b/city_metrix/layers/average_net_building_height.py index 2a215c7..d0f49f2 100644 --- a/city_metrix/layers/average_net_building_height.py +++ b/city_metrix/layers/average_net_building_height.py @@ -6,9 +6,14 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection class AverageNetBuildingHeight(Layer): - def __init__(self, scale_meters=100, **kwargs): + """ + Attributes: + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, spatial_resolution=100, **kwargs): super().__init__(**kwargs) - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): # https://ghsl.jrc.ec.europa.eu/ghs_buH2023.php @@ -19,7 +24,7 @@ def get_data(self, bbox): anbh = ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A") data = (get_image_collection( - ee.ImageCollection(anbh), bbox, self.scale_meters, "average net building height") + ee.ImageCollection(anbh), bbox, self.spatial_resolution, "average net building height") .b1) return data diff --git a/city_metrix/layers/built_up_height.py b/city_metrix/layers/built_up_height.py index bc89d70..ab080f5 100644 --- a/city_metrix/layers/built_up_height.py +++ b/city_metrix/layers/built_up_height.py @@ -7,9 +7,14 @@ class BuiltUpHeight(Layer): - def __init__(self, scale_meters=100, **kwargs): + """ + Attributes: + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, spatial_resolution=100, **kwargs): super().__init__(**kwargs) - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): # Notes for Heat project: @@ -19,6 +24,6 @@ def get_data(self, bbox): # ee.ImageCollection("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_R2023A") built_height = ee.Image("JRC/GHSL/P2023A/GHS_BUILT_H/2018") - data = get_image_collection(ee.ImageCollection(built_height), bbox, self.scale_meters, "built up height") + data = get_image_collection(ee.ImageCollection(built_height), bbox, self.spatial_resolution, "built up height") return data.built_height diff --git a/city_metrix/layers/esa_world_cover.py b/city_metrix/layers/esa_world_cover.py index b306d8c..a4b0f65 100644 --- a/city_metrix/layers/esa_world_cover.py +++ b/city_metrix/layers/esa_world_cover.py @@ -20,29 +20,36 @@ class EsaWorldCoverClass(Enum): MOSS_AND_LICHEN = 100 class EsaWorldCover(Layer): + """ + Attributes: + land_cover_class: Enum value from EsaWorldCoverClass + year: year used for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + STAC_CATALOG_URI = "https://services.terrascope.be/stac/" STAC_COLLECTION_ID = "urn:eop:VITO:ESA_WorldCover_10m_2020_AWS_V1" STAC_ASSET_ID = "ESA_WORLDCOVER_10M_MAP" - def __init__(self, land_cover_class=None, year=2020, scale_meters=10, **kwargs): + def __init__(self, land_cover_class=None, year=2020, spatial_resolution=10, **kwargs): super().__init__(**kwargs) self.land_cover_class = land_cover_class self.year = year - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): if self.year == 2020: data = get_image_collection( ee.ImageCollection("ESA/WorldCover/v100"), bbox, - self.scale_meters, + self.spatial_resolution, "ESA world cover" ).Map elif self.year == 2021: data = get_image_collection( ee.ImageCollection("ESA/WorldCover/v200"), bbox, - self.scale_meters, + self.spatial_resolution, "ESA world cover" ).Map diff --git a/city_metrix/layers/land_surface_temperature.py b/city_metrix/layers/land_surface_temperature.py index da4922c..931cb2e 100644 --- a/city_metrix/layers/land_surface_temperature.py +++ b/city_metrix/layers/land_surface_temperature.py @@ -6,11 +6,18 @@ import xarray class LandSurfaceTemperature(Layer): - def __init__(self, start_date="2013-01-01", end_date="2023-01-01", scale_meters=30, **kwargs): + """ + Attributes: + start_date: starting date for data retrieval + end_date: ending date for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resolution=30, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): def cloud_mask(image): @@ -31,5 +38,5 @@ def apply_scale_factors(image): .map(apply_scale_factors) \ .reduce(ee.Reducer.mean()) - data = get_image_collection(ee.ImageCollection(l8_st), bbox, self.scale_meters, "LST").ST_B10_mean + data = get_image_collection(ee.ImageCollection(l8_st), bbox, self.spatial_resolution, "LST").ST_B10_mean return data diff --git a/city_metrix/layers/nasa_dem.py b/city_metrix/layers/nasa_dem.py index 35daa0c..b5ac45d 100644 --- a/city_metrix/layers/nasa_dem.py +++ b/city_metrix/layers/nasa_dem.py @@ -6,9 +6,14 @@ class NasaDEM(Layer): - def __init__(self, scale_meters=30, **kwargs): + """ + Attributes: + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, spatial_resolution=30, **kwargs): super().__init__(**kwargs) - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): dataset = ee.Image("NASA/NASADEM_HGT/001") @@ -17,6 +22,6 @@ def get_data(self, bbox): .select('elevation') .mean() ) - data = get_image_collection(nasa_dem, bbox, self.scale_meters, "NASA DEM").elevation + data = get_image_collection(nasa_dem, bbox, self.spatial_resolution, "NASA DEM").elevation return data diff --git a/city_metrix/layers/natural_areas.py b/city_metrix/layers/natural_areas.py index 61c648a..5efe4a8 100644 --- a/city_metrix/layers/natural_areas.py +++ b/city_metrix/layers/natural_areas.py @@ -5,12 +5,17 @@ from .esa_world_cover import EsaWorldCover, EsaWorldCoverClass class NaturalAreas(Layer): - def __init__(self, scale_meters=10, **kwargs): + """ + Attributes: + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, spatial_resolution=10, **kwargs): super().__init__(**kwargs) - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): - esa_world_cover = EsaWorldCover(scale_meters=self.scale_meters).get_data(bbox) + esa_world_cover = EsaWorldCover(spatial_resolution=self.spatial_resolution).get_data(bbox) reclass_map = { EsaWorldCoverClass.TREE_COVER.value: 1, EsaWorldCoverClass.SHRUBLAND.value: 1, diff --git a/city_metrix/layers/ndvi_sentinel2_gee.py b/city_metrix/layers/ndvi_sentinel2_gee.py index ca7fc05..89b4c09 100644 --- a/city_metrix/layers/ndvi_sentinel2_gee.py +++ b/city_metrix/layers/ndvi_sentinel2_gee.py @@ -4,16 +4,18 @@ class NdviSentinel2(Layer): """" NDVI = Sentinel-2 Normalized Difference Vegetation Index - param: year: The satellite imaging year. + Attributes: + year: The satellite imaging year. + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) return: a rioxarray-format DataArray Author of associated Jupyter notebook: Ted.Wong@wri.org Notebook: https://github.com/wri/cities-cities4forests-indicators/blob/dev-eric/scripts/extract-VegetationCover.ipynb Reference: https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index """ - def __init__(self, year=None, scale_meters=10, **kwargs): + def __init__(self, year=None, spatial_resolution=10, **kwargs): super().__init__(**kwargs) self.year = year - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): if self.year is None: @@ -40,7 +42,7 @@ def calculate_ndvi(image): ndvi_mosaic = ndvi.qualityMosaic('NDVI') ic = ee.ImageCollection(ndvi_mosaic) - ndvi_data = get_image_collection(ic, bbox, self.scale_meters, "NDVI") + ndvi_data = get_image_collection(ic, bbox, self.spatial_resolution, "NDVI") xdata = ndvi_data.to_dataarray() diff --git a/city_metrix/layers/tree_canopy_height.py b/city_metrix/layers/tree_canopy_height.py index 9141f84..fee1697 100644 --- a/city_metrix/layers/tree_canopy_height.py +++ b/city_metrix/layers/tree_canopy_height.py @@ -6,12 +6,17 @@ import ee class TreeCanopyHeight(Layer): + """ + Attributes: + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + name = "tree_canopy_height" NO_DATA_VALUE = 0 - def __init__(self, scale_meters=1, **kwargs): + def __init__(self, spatial_resolution=1, **kwargs): super().__init__(**kwargs) - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): canopy_ht = ee.ImageCollection("projects/meta-forest-monitoring-okw37/assets/CanopyHeight") @@ -19,6 +24,6 @@ def get_data(self, bbox): canopy_ht = canopy_ht.reduce(ee.Reducer.mean()).rename("cover_code") data = get_image_collection(ee.ImageCollection(canopy_ht), bbox, - self.scale_meters, "tree canopy height") + self.spatial_resolution, "tree canopy height") return data.cover_code diff --git a/city_metrix/layers/tree_cover.py b/city_metrix/layers/tree_cover.py index 9caa483..98bc481 100644 --- a/city_metrix/layers/tree_cover.py +++ b/city_metrix/layers/tree_cover.py @@ -8,15 +8,19 @@ class TreeCover(Layer): """ Merged tropical and nontropical tree cover from WRI + Attributes: + min_tree_cover: minimum tree-cover values used for filtering results + max_tree_cover: maximum tree-cover values used for filtering results + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) """ - + NO_DATA_VALUE = 255 - def __init__(self, min_tree_cover=None, max_tree_cover=None, scale_meters=10, **kwargs): + def __init__(self, min_tree_cover=None, max_tree_cover=None, spatial_resolution=10, **kwargs): super().__init__(**kwargs) self.min_tree_cover = min_tree_cover self.max_tree_cover = max_tree_cover - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): tropics = ee.ImageCollection('projects/wri-datalab/TropicalTreeCover') @@ -24,7 +28,7 @@ def get_data(self, bbox): merged_ttc = tropics.merge(nontropics) ttc_image = merged_ttc.reduce(ee.Reducer.mean()).rename('ttc') - data = get_image_collection(ee.ImageCollection(ttc_image), bbox, self.scale_meters, "tree cover").ttc + data = get_image_collection(ee.ImageCollection(ttc_image), bbox, self.spatial_resolution, "tree cover").ttc if self.min_tree_cover is not None: data = data.where(data >= self.min_tree_cover) diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py index da83219..fe69c75 100644 --- a/city_metrix/layers/urban_land_use.py +++ b/city_metrix/layers/urban_land_use.py @@ -7,10 +7,16 @@ class UrbanLandUse(Layer): - def __init__(self, band='lulc', scale_meters=5, **kwargs): + """ + Attributes: + band: raster band used for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, band='lulc', spatial_resolution=5, **kwargs): super().__init__(**kwargs) self.band = band - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): dataset = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1") @@ -28,6 +34,6 @@ def get_data(self, bbox): .rename('lulc') ) - data = get_image_collection(ulu, bbox, self.scale_meters, "urban land use").lulc + data = get_image_collection(ulu, bbox, self.spatial_resolution, "urban land use").lulc return data diff --git a/city_metrix/layers/world_pop.py b/city_metrix/layers/world_pop.py index d77cc2e..dff494a 100644 --- a/city_metrix/layers/world_pop.py +++ b/city_metrix/layers/world_pop.py @@ -6,9 +6,14 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection class WorldPop(Layer): - def __init__(self, scale_meters=100, **kwargs): + """ + Attributes: + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ + + def __init__(self, spatial_resolution=100, **kwargs): super().__init__(**kwargs) - self.scale_meters = scale_meters + self.spatial_resolution = spatial_resolution def get_data(self, bbox): # load population @@ -20,5 +25,5 @@ def get_data(self, bbox): .mean() ) - data = get_image_collection(world_pop, bbox, self.scale_meters, "world pop") + data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop") return data.population diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index bef67ec..ce220a5 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,93 +1,106 @@ from city_metrix.layers import ( + Layer, Albedo, AlosDSM, AverageNetBuildingHeight, + BuiltUpHeight, EsaWorldCover, EsaWorldCoverClass, LandSurfaceTemperature, NasaDEM, NaturalAreas, + NdviSentinel2, TreeCanopyHeight, TreeCover, UrbanLandUse, - WorldPop, NdviSentinel2, BuiltUpHeight + WorldPop ) from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 +""" +Note: To add a test for another scalable layer that has the spatial_resolution property: +1. Add the class name to the city_metrix.layers import statement above +2. Specify a minimal class instance in the set below. Do no specify the spatial_resolution + property in the instance definition. +""" +CLASSES_WITH_spatial_resolution_PROPERTY = \ + { + 'Albedo()', + 'AlosDSM()', + 'AverageNetBuildingHeight()', + 'BuiltUpHeight()', + 'EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)', + 'LandSurfaceTemperature()', + 'NasaDEM()', + 'NaturalAreas()', + 'NdviSentinel2(year=2023)', + 'TreeCanopyHeight()', + 'TreeCover()', + 'UrbanLandUse()', + 'WorldPop()' + } + COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 -def test_albedo_scale(): - doubled_default_resolution = 2 * Albedo().scale_meters - layer = Albedo(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_alos_dsm_scale(): - doubled_default_resolution = 2 * AlosDSM().scale_meters - layer = AlosDSM(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_average_net_building_height_scale(): - doubled_default_resolution = 2 * AverageNetBuildingHeight().scale_meters - layer = AverageNetBuildingHeight(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_built_up_height_scale(): - doubled_default_resolution = 2 * BuiltUpHeight().scale_meters - layer = BuiltUpHeight(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_esa_world_cover_scale(): - doubled_default_resolution = 2 * EsaWorldCover().scale_meters - layer = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP, scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -# TODO -# def test_high_land_surface_temperature_scale(): - -def test_land_surface_temperature_scale(): - doubled_default_resolution = 2 * LandSurfaceTemperature().scale_meters - layer = LandSurfaceTemperature(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_nasa_dem_scale(): - doubled_default_resolution = 2 * NasaDEM().scale_meters - layer = NasaDEM(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_natural_areas_scale(): - doubled_default_resolution = 2 * NaturalAreas().scale_meters - layer = NaturalAreas(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_ndvi_sentinel2_scale(): - doubled_default_resolution = 2 * NdviSentinel2().scale_meters - layer = NdviSentinel2(year=2023, scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -# TODO -# def test_smart_surface_lulc_scale(): - -def test_tree_canopy_height(): - doubled_default_resolution = 2 * TreeCanopyHeight().scale_meters - layer = TreeCanopyHeight(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_tree_cover(): - doubled_default_resolution = 2 * TreeCover().scale_meters - layer = TreeCover(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_urban_land_use(): - doubled_default_resolution = 2 * UrbanLandUse().scale_meters - layer = UrbanLandUse(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - -def test_world_pop(): - doubled_default_resolution = 2 * WorldPop().scale_meters - layer = WorldPop(scale_meters=doubled_default_resolution) - evaluate_layer(layer, doubled_default_resolution) - +def test_scale__meters_property_for_all_scalable_layers(): + for class_instance_str in CLASSES_WITH_spatial_resolution_PROPERTY: + is_valid, except_str = validate_layer_instance(class_instance_str) + if is_valid is False: + raise Exception(except_str) + + class_instance = eval(class_instance_str) + + # Double the spatial_resolution for the specified Class + doubled_default_resolution = 2 * class_instance.spatial_resolution + class_instance.spatial_resolution=doubled_default_resolution + + evaluate_layer(class_instance, doubled_default_resolution) + + +def test_function_validate_layer_instance(): + is_valid, except_str = validate_layer_instance(Albedo()) + assert is_valid is False + is_valid, except_str = validate_layer_instance('t') + assert is_valid is False + is_valid, except_str = validate_layer_instance('Layer()') + assert is_valid is False + is_valid, except_str = validate_layer_instance('OpenStreetMap()') + assert is_valid is False + is_valid, except_str = validate_layer_instance('Albedo(spatial_resolution = 2)') + assert is_valid is False + +def validate_layer_instance(obj_string): + is_valid = True + except_str = None + obj_eval = None + + if not type(obj_string) == str: + is_valid = False + except_str = "Specified object '%s' must be specified as a string." % obj_string + return is_valid, except_str + + try: + obj_eval = eval(obj_string) + except: + is_valid = False + except_str = "Specified object '%s' is not a class instance." % obj_string + return is_valid, except_str + + if not type(obj_eval).__bases__[0] == Layer: + is_valid = False + except_str = "Specified object '%s' is not a valid Layer class instance." % obj_string + elif not hasattr(obj_eval, 'spatial_resolution'): + is_valid = False + except_str = "Specified class '%s' does not have the spatial_resolution property." % obj_string + elif not obj_string.find('spatial_resolution') == -1: + is_valid = False + except_str = "Do not specify spatial_resolution property value in object '%s'." % obj_string + elif obj_eval.spatial_resolution is None: + is_valid = False + except_str = "Class signature cannot specify None for default value for class." + + return is_valid, except_str def evaluate_layer(layer, expected_resolution): data = layer.get_data(BBOX) From c33c58cb3f609b68601aaa82a21e32891617efe1 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 29 Aug 2024 15:37:39 -0700 Subject: [PATCH 03/32] renamed function and trying to resolve GH Action failure. --- tests/test_layer_parameters.py | 35 +++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index ce220a5..305523b 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -25,25 +25,25 @@ """ CLASSES_WITH_spatial_resolution_PROPERTY = \ { - 'Albedo()', - 'AlosDSM()', - 'AverageNetBuildingHeight()', - 'BuiltUpHeight()', + # 'Albedo()', + # 'AlosDSM()', + # 'AverageNetBuildingHeight()', + # 'BuiltUpHeight()', 'EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)', - 'LandSurfaceTemperature()', - 'NasaDEM()', - 'NaturalAreas()', - 'NdviSentinel2(year=2023)', - 'TreeCanopyHeight()', - 'TreeCover()', - 'UrbanLandUse()', - 'WorldPop()' + # 'LandSurfaceTemperature()', + # 'NasaDEM()', + # 'NaturalAreas()', + # 'NdviSentinel2(year=2023)', + # 'TreeCanopyHeight()', + # 'TreeCover()', + # 'UrbanLandUse()', + # 'WorldPop()' } COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 -def test_scale__meters_property_for_all_scalable_layers(): +def test_spatial_resolution_for_all_scalable_layers(): for class_instance_str in CLASSES_WITH_spatial_resolution_PROPERTY: is_valid, except_str = validate_layer_instance(class_instance_str) if is_valid is False: @@ -104,13 +104,14 @@ def validate_layer_instance(obj_string): def evaluate_layer(layer, expected_resolution): data = layer.get_data(BBOX) - est_actual_resolution = get_resolution_estimate(data) - assert expected_resolution == est_actual_resolution + actual_estimated_resolution, y_cells, y_diff = get_resolution_estimate(data) + print ('y_cells %s, y_diff %s' % (y_cells, y_diff)) + assert expected_resolution == actual_estimated_resolution def get_resolution_estimate(data): - y_cells = data['y'].size - 1 + y_cells = float(data['y'].size - 1) y_min = data['y'].values.min() y_max = data['y'].values.max() y_diff = y_max - y_min ry = round(y_diff/y_cells) - return ry + return ry, y_cells, y_diff From c476aa110ac926b38ac31bd4e1d5bae0ac20a2f1 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 29 Aug 2024 15:51:58 -0700 Subject: [PATCH 04/32] Exposing parameters to understand Action failure --- tests/test_layer_parameters.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 305523b..7ded19b 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -104,8 +104,8 @@ def validate_layer_instance(obj_string): def evaluate_layer(layer, expected_resolution): data = layer.get_data(BBOX) - actual_estimated_resolution, y_cells, y_diff = get_resolution_estimate(data) - print ('y_cells %s, y_diff %s' % (y_cells, y_diff)) + actual_estimated_resolution, y_cells, y_min, y_max = get_resolution_estimate(data) + print ('y_cells %s, y_min %s, y_max %s' % (y_cells, y_min, y_max)) assert expected_resolution == actual_estimated_resolution def get_resolution_estimate(data): @@ -114,4 +114,4 @@ def get_resolution_estimate(data): y_max = data['y'].values.max() y_diff = y_max - y_min ry = round(y_diff/y_cells) - return ry, y_cells, y_diff + return ry, y_cells, y_min, y_max From b532b25bbee176c1817843ebafaec0b80208a8e0 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 29 Aug 2024 18:53:39 -0700 Subject: [PATCH 05/32] Handle other coordinate units --- .../conftest.py | 2 +- tests/test_layer_parameters.py | 51 ++++++++++++------- tests/tools/__init__.py | 0 {tools => tests/tools}/general_tools.py | 0 tests/tools/spatial_tools.py | 16 ++++++ 5 files changed, 50 insertions(+), 19 deletions(-) create mode 100644 tests/tools/__init__.py rename {tools => tests/tools}/general_tools.py (100%) create mode 100644 tests/tools/spatial_tools.py diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py index 5882053..f8727f5 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py +++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py @@ -5,7 +5,7 @@ from collections import namedtuple from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 -from tools.general_tools import create_target_folder, is_valid_path +from tests.tools.general_tools import create_target_folder, is_valid_path # RUN_DUMPS is the master control for whether the writes and tests are executed # Setting RUN_DUMPS to True turns on code execution. diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 7ded19b..1b91200 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,21 +1,23 @@ +from pyproj import CRS from city_metrix.layers import ( Layer, - Albedo, - AlosDSM, - AverageNetBuildingHeight, - BuiltUpHeight, + # Albedo, + # AlosDSM, + # AverageNetBuildingHeight, + # BuiltUpHeight, EsaWorldCover, EsaWorldCoverClass, - LandSurfaceTemperature, - NasaDEM, - NaturalAreas, - NdviSentinel2, - TreeCanopyHeight, - TreeCover, - UrbanLandUse, - WorldPop + # LandSurfaceTemperature, + # NasaDEM, + # NaturalAreas, + # NdviSentinel2, + # TreeCanopyHeight, + # TreeCover, + # UrbanLandUse, + # WorldPop ) from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 +from tests.tools.spatial_tools import get_distance_between_geocoordinates """ Note: To add a test for another scalable layer that has the spatial_resolution property: @@ -104,14 +106,27 @@ def validate_layer_instance(obj_string): def evaluate_layer(layer, expected_resolution): data = layer.get_data(BBOX) - actual_estimated_resolution, y_cells, y_min, y_max = get_resolution_estimate(data) - print ('y_cells %s, y_min %s, y_max %s' % (y_cells, y_min, y_max)) + actual_estimated_resolution, uu = get_spatial_resolution_estimate(data) + print ('uu %s' % uu) assert expected_resolution == actual_estimated_resolution -def get_resolution_estimate(data): +def get_spatial_resolution_estimate(data): y_cells = float(data['y'].size - 1) + y_min = data['y'].values.min() y_max = data['y'].values.max() - y_diff = y_max - y_min - ry = round(y_diff/y_cells) - return ry, y_cells, y_min, y_max + + crs = CRS.from_string(data.crs) + uu = crs.axis_info[0].unit_name + if uu == 'metre': + y_diff = y_max - y_min + else: + lat1 = y_min + lat2 = y_min + lon1 = data['x'].values.min() + lon2 = data['x'].values.max() + y_diff = get_distance_between_geocoordinates(lat1, lon1, lat2, lon2) + + ry = round(y_diff / y_cells) + + return ry, uu diff --git a/tests/tools/__init__.py b/tests/tools/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tools/general_tools.py b/tests/tools/general_tools.py similarity index 100% rename from tools/general_tools.py rename to tests/tools/general_tools.py diff --git a/tests/tools/spatial_tools.py b/tests/tools/spatial_tools.py new file mode 100644 index 0000000..73ab093 --- /dev/null +++ b/tests/tools/spatial_tools.py @@ -0,0 +1,16 @@ +import math + +EARTH_RADIUS = 6378.137 # Radius of earth in KM + +def get_distance_between_geocoordinates(lat1, lon1, lat2, lon2): + dLat = lat2 * math.pi / 180 - lat1 * math.pi / 180 + dLon = lon2 * math.pi / 180 - lon1 * math.pi / 180 + a = math.sin(dLat/2) * math.sin(dLat/2) + \ + math.cos(lat1 * math.pi / 180) * \ + math.cos(lat2 * math.pi / 180) * \ + math.sin(dLon/2) * math.sin(dLon/2) + c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) + d = EARTH_RADIUS * c + return d * 1000 # meters + + From fecab9766ca66eed74a94ac19f73d74e0aac0bc5 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 29 Aug 2024 19:15:04 -0700 Subject: [PATCH 06/32] Fixed error in location of test tools --- tests/test_layer_dimensions.py | 2 +- tests/test_layer_parameters.py | 24 ++++++++++----------- tests/tools.py | 35 ------------------------------- tests/tools/general_tools.py | 38 ++++++++++++++++++++++++++++++++-- 4 files changed, 49 insertions(+), 50 deletions(-) delete mode 100644 tests/tools.py diff --git a/tests/test_layer_dimensions.py b/tests/test_layer_dimensions.py index ee32af2..c3fdd6d 100644 --- a/tests/test_layer_dimensions.py +++ b/tests/test_layer_dimensions.py @@ -1,6 +1,6 @@ from city_metrix.layers import NdviSentinel2 from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 -from tests.tools import post_process_layer +from tests.tools.general_tools import post_process_layer COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 1b91200..57d7351 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,20 +1,20 @@ from pyproj import CRS from city_metrix.layers import ( Layer, - # Albedo, - # AlosDSM, - # AverageNetBuildingHeight, - # BuiltUpHeight, + Albedo, + AlosDSM, + AverageNetBuildingHeight, + BuiltUpHeight, EsaWorldCover, EsaWorldCoverClass, - # LandSurfaceTemperature, - # NasaDEM, - # NaturalAreas, - # NdviSentinel2, - # TreeCanopyHeight, - # TreeCover, - # UrbanLandUse, - # WorldPop + LandSurfaceTemperature, + NasaDEM, + NaturalAreas, + NdviSentinel2, + TreeCanopyHeight, + TreeCover, + UrbanLandUse, + WorldPop ) from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 from tests.tools.spatial_tools import get_distance_between_geocoordinates diff --git a/tests/tools.py b/tests/tools.py deleted file mode 100644 index 99425df..0000000 --- a/tests/tools.py +++ /dev/null @@ -1,35 +0,0 @@ -import numpy as np - -def post_process_layer(data, value_threshold=0.4, convert_to_percentage=True): - """ - Applies the standard post-processing adjustment used for rendering of NDVI including masking - to a threshold and conversion to percentage values. - :param value_threshold: (float) minimum threshold for keeping values - :param convert_to_percentage: (bool) controls whether NDVI values are converted to a percentage - :return: A rioxarray-format DataArray - """ - # Remove values less than the specified threshold - if value_threshold is not None: - data = data.where(data >= value_threshold) - - # Convert to percentage in byte data_type - if convert_to_percentage is True: - data = convert_ratio_to_percentage(data) - - return data - -def convert_ratio_to_percentage(data): - """ - Converts xarray variable from a ratio to a percentage - :param data: (xarray) xarray to be converted - :return: A rioxarray-format DataArray - """ - - # convert to percentage and to bytes for efficient storage - values_as_percent = np.round(data * 100).astype(np.uint8) - - # reset CRS - source_crs = data.rio.crs - values_as_percent.rio.write_crs(source_crs, inplace=True) - - return values_as_percent diff --git a/tests/tools/general_tools.py b/tests/tools/general_tools.py index b38d1b5..99bd894 100644 --- a/tests/tools/general_tools.py +++ b/tests/tools/general_tools.py @@ -1,6 +1,6 @@ import os import tempfile - +import numpy as np def is_valid_path(path: str): return os.path.exists(path) @@ -20,4 +20,38 @@ def remove_all_files_in_directory(directory): if os.path.isfile(file_path): os.remove(file_path) except Exception as e: - print(f"Error: {e}") \ No newline at end of file + print(f"Error: {e}") + +def post_process_layer(data, value_threshold=0.4, convert_to_percentage=True): + """ + Applies the standard post-processing adjustment used for rendering of NDVI including masking + to a threshold and conversion to percentage values. + :param value_threshold: (float) minimum threshold for keeping values + :param convert_to_percentage: (bool) controls whether NDVI values are converted to a percentage + :return: A rioxarray-format DataArray + """ + # Remove values less than the specified threshold + if value_threshold is not None: + data = data.where(data >= value_threshold) + + # Convert to percentage in byte data_type + if convert_to_percentage is True: + data = convert_ratio_to_percentage(data) + + return data + +def convert_ratio_to_percentage(data): + """ + Converts xarray variable from a ratio to a percentage + :param data: (xarray) xarray to be converted + :return: A rioxarray-format DataArray + """ + + # convert to percentage and to bytes for efficient storage + values_as_percent = np.round(data * 100).astype(np.uint8) + + # reset CRS + source_crs = data.rio.crs + values_as_percent.rio.write_crs(source_crs, inplace=True) + + return values_as_percent \ No newline at end of file From c820296731c6a77c885d39f023389aa8cae60d89 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 29 Aug 2024 19:50:33 -0700 Subject: [PATCH 07/32] Changed estimate of spatial resolution --- tests/test_layer_parameters.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 57d7351..043d917 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -106,8 +106,7 @@ def validate_layer_instance(obj_string): def evaluate_layer(layer, expected_resolution): data = layer.get_data(BBOX) - actual_estimated_resolution, uu = get_spatial_resolution_estimate(data) - print ('uu %s' % uu) + actual_estimated_resolution = get_spatial_resolution_estimate(data) assert expected_resolution == actual_estimated_resolution def get_spatial_resolution_estimate(data): @@ -117,16 +116,21 @@ def get_spatial_resolution_estimate(data): y_max = data['y'].values.max() crs = CRS.from_string(data.crs) - uu = crs.axis_info[0].unit_name - if uu == 'metre': + crs_units = crs.axis_info[0].unit_name + if crs_units == 'metre': y_diff = y_max - y_min - else: + elif crs_units == 'foot': + feet_to_meter = 0.3048 + y_diff = (y_max - y_min) * feet_to_meter + elif crs_units == 'degree': lat1 = y_min - lat2 = y_min + lat2 = y_max lon1 = data['x'].values.min() - lon2 = data['x'].values.max() + lon2 = lon1 y_diff = get_distance_between_geocoordinates(lat1, lon1, lat2, lon2) + else: + raise Exception('Unhandled projection units: %s' % crs_units) ry = round(y_diff / y_cells) - return ry, uu + return ry From ea8cf8855e4f53cecdf6329ea3a63765ae0131e7 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 29 Aug 2024 22:14:42 -0700 Subject: [PATCH 08/32] Changed object definition --- tests/test_layer_parameters.py | 110 ++++++++++++++++----------------- 1 file changed, 53 insertions(+), 57 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 043d917..430703f 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -5,8 +5,7 @@ AlosDSM, AverageNetBuildingHeight, BuiltUpHeight, - EsaWorldCover, - EsaWorldCoverClass, + EsaWorldCover, EsaWorldCoverClass, LandSurfaceTemperature, NasaDEM, NaturalAreas, @@ -14,7 +13,7 @@ TreeCanopyHeight, TreeCover, UrbanLandUse, - WorldPop + WorldPop, OpenStreetMap ) from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 from tests.tools.spatial_tools import get_distance_between_geocoordinates @@ -25,98 +24,95 @@ 2. Specify a minimal class instance in the set below. Do no specify the spatial_resolution property in the instance definition. """ -CLASSES_WITH_spatial_resolution_PROPERTY = \ +CLASSES_WITH_SPATIAL_RESOLUTION_PROPERTY = \ { - # 'Albedo()', - # 'AlosDSM()', - # 'AverageNetBuildingHeight()', - # 'BuiltUpHeight()', - 'EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)', - # 'LandSurfaceTemperature()', - # 'NasaDEM()', - # 'NaturalAreas()', - # 'NdviSentinel2(year=2023)', - # 'TreeCanopyHeight()', - # 'TreeCover()', - # 'UrbanLandUse()', - # 'WorldPop()' + Albedo(), + AlosDSM(), + AverageNetBuildingHeight(), + BuiltUpHeight(), + EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP), + LandSurfaceTemperature(), + NasaDEM(), + NaturalAreas(), + NdviSentinel2(year=2023), + TreeCanopyHeight(), + TreeCover(), + UrbanLandUse(), + WorldPop() } COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 def test_spatial_resolution_for_all_scalable_layers(): - for class_instance_str in CLASSES_WITH_spatial_resolution_PROPERTY: - is_valid, except_str = validate_layer_instance(class_instance_str) + for obj in CLASSES_WITH_SPATIAL_RESOLUTION_PROPERTY: + is_valid, except_str = validate_layer_instance(obj) if is_valid is False: raise Exception(except_str) - class_instance = eval(class_instance_str) + cls = get_class_from_instance(obj) # Double the spatial_resolution for the specified Class - doubled_default_resolution = 2 * class_instance.spatial_resolution - class_instance.spatial_resolution=doubled_default_resolution + doubled_default_resolution = 2 * cls.spatial_resolution + obj.spatial_resolution=doubled_default_resolution - evaluate_layer(class_instance, doubled_default_resolution) + evaluate_layer(obj, doubled_default_resolution) def test_function_validate_layer_instance(): - is_valid, except_str = validate_layer_instance(Albedo()) - assert is_valid is False is_valid, except_str = validate_layer_instance('t') assert is_valid is False - is_valid, except_str = validate_layer_instance('Layer()') + is_valid, except_str = validate_layer_instance(EsaWorldCoverClass.BUILT_UP) assert is_valid is False - is_valid, except_str = validate_layer_instance('OpenStreetMap()') + is_valid, except_str = validate_layer_instance(OpenStreetMap()) assert is_valid is False - is_valid, except_str = validate_layer_instance('Albedo(spatial_resolution = 2)') + is_valid, except_str = validate_layer_instance(Albedo(spatial_resolution = 2)) assert is_valid is False -def validate_layer_instance(obj_string): +def validate_layer_instance(obj): is_valid = True except_str = None - obj_eval = None - - if not type(obj_string) == str: - is_valid = False - except_str = "Specified object '%s' must be specified as a string." % obj_string - return is_valid, except_str - - try: - obj_eval = eval(obj_string) - except: - is_valid = False - except_str = "Specified object '%s' is not a class instance." % obj_string - return is_valid, except_str - if not type(obj_eval).__bases__[0] == Layer: - is_valid = False - except_str = "Specified object '%s' is not a valid Layer class instance." % obj_string - elif not hasattr(obj_eval, 'spatial_resolution'): + if not obj.__class__.__bases__[0] == Layer: is_valid = False - except_str = "Specified class '%s' does not have the spatial_resolution property." % obj_string - elif not obj_string.find('spatial_resolution') == -1: - is_valid = False - except_str = "Do not specify spatial_resolution property value in object '%s'." % obj_string - elif obj_eval.spatial_resolution is None: - is_valid = False - except_str = "Class signature cannot specify None for default value for class." + except_str = "Specified object '%s' is not a valid Layer class instance." % obj + else: + cls = get_class_from_instance(obj) + cls_name = type(cls).__name__ + if not hasattr(obj, 'spatial_resolution'): + is_valid = False + except_str = "Class '%s' does not have spatial_resolution property." % cls_name + elif not obj.spatial_resolution == cls.spatial_resolution: + is_valid = False + except_str = "Do not specify spatial_resolution property value for class '%s'." % cls_name + elif cls.spatial_resolution is None: + is_valid = False + except_str = "Signature of class %s must specify a non-null default value for spatial_resolution. Please correct." % cls_name return is_valid, except_str +def get_class_from_instance(obj): + cls = obj.__class__() + return cls + def evaluate_layer(layer, expected_resolution): data = layer.get_data(BBOX) - actual_estimated_resolution = get_spatial_resolution_estimate(data) + actual_estimated_resolution = estimate_spatial_resolution(data) assert expected_resolution == actual_estimated_resolution -def get_spatial_resolution_estimate(data): +def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) y_min = data['y'].values.min() y_max = data['y'].values.max() - crs = CRS.from_string(data.crs) - crs_units = crs.axis_info[0].unit_name + try: + crs = CRS.from_string(data.crs) + crs_units = crs.axis_info[0].unit_name + except: + # if xarray doesn't have crs property, assume meters + crs_units = 'metre' + if crs_units == 'metre': y_diff = y_max - y_min elif crs_units == 'foot': From 08e3d8b25c2113635cfd15fb1a7603a462a26f65 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 12:53:27 -0700 Subject: [PATCH 09/32] Restructured tests to have individual test per layer --- tests/test_layer_parameters.py | 124 +++++++++++++++++++++++---------- 1 file changed, 87 insertions(+), 37 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 430703f..c939202 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -19,44 +19,99 @@ from tests.tools.spatial_tools import get_distance_between_geocoordinates """ -Note: To add a test for another scalable layer that has the spatial_resolution property: +Evaluation of spatial_resolution property +To add a test for a scalable layer that has the spatial_resolution property: 1. Add the class name to the city_metrix.layers import statement above -2. Specify a minimal class instance in the set below. Do no specify the spatial_resolution - property in the instance definition. +2. Copy an existing test_*_spatial_resolution() test + a. rename for the new layer + b. specify a minimal class instance for the layer, not specifying the spatial_resolution attribute. """ -CLASSES_WITH_SPATIAL_RESOLUTION_PROPERTY = \ - { - Albedo(), - AlosDSM(), - AverageNetBuildingHeight(), - BuiltUpHeight(), - EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP), - LandSurfaceTemperature(), - NasaDEM(), - NaturalAreas(), - NdviSentinel2(year=2023), - TreeCanopyHeight(), - TreeCover(), - UrbanLandUse(), - WorldPop() - } COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 -def test_spatial_resolution_for_all_scalable_layers(): - for obj in CLASSES_WITH_SPATIAL_RESOLUTION_PROPERTY: - is_valid, except_str = validate_layer_instance(obj) - if is_valid is False: - raise Exception(except_str) - - cls = get_class_from_instance(obj) - - # Double the spatial_resolution for the specified Class - doubled_default_resolution = 2 * cls.spatial_resolution - obj.spatial_resolution=doubled_default_resolution +def test_albedo_spatial_resolution(): + class_instance = Albedo() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_alos_dsm_spatial_resolution(): + class_instance = AlosDSM() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_average_net_building_height_spatial_resolution(): + class_instance = AverageNetBuildingHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_built_up_height_spatial_resolution(): + class_instance = BuiltUpHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_esa_world_cover_spatial_resolution(): + class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_land_surface_temperature_spatial_resolution(): + class_instance = LandSurfaceTemperature() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_nasa_dem_spatial_resolution(): + class_instance = NasaDEM() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_natural_areas_spatial_resolution(): + class_instance = NaturalAreas() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_ndvi_sentinel2_spatial_resolution(): + class_instance = NdviSentinel2(year=2023) + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_tree_canopy_height_spatial_resolution(): + class_instance = TreeCanopyHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_tree_cover_spatial_resolution(): + class_instance = TreeCover() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_urban_land_use_spatial_resolution(): + class_instance = UrbanLandUse() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_world_pop_spatial_resolution(): + class_instance = WorldPop() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def evaluate_resolution__property(obj): + is_valid, except_str = validate_layer_instance(obj) + if is_valid is False: + raise Exception(except_str) + + cls = get_class_from_instance(obj) + + # Double the spatial_resolution for the specified Class + doubled_default_resolution = 2 * cls.spatial_resolution + obj.spatial_resolution=doubled_default_resolution + + data = obj.get_data(BBOX) + + expected_resolution = doubled_default_resolution + actual_estimated_resolution = estimate_spatial_resolution(data) - evaluate_layer(obj, doubled_default_resolution) + return expected_resolution, actual_estimated_resolution def test_function_validate_layer_instance(): @@ -95,11 +150,6 @@ def get_class_from_instance(obj): cls = obj.__class__() return cls -def evaluate_layer(layer, expected_resolution): - data = layer.get_data(BBOX) - actual_estimated_resolution = estimate_spatial_resolution(data) - assert expected_resolution == actual_estimated_resolution - def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) @@ -113,7 +163,7 @@ def estimate_spatial_resolution(data): # if xarray doesn't have crs property, assume meters crs_units = 'metre' - if crs_units == 'metre': + if crs_units == 'metre' or crs_units == 'm': y_diff = y_max - y_min elif crs_units == 'foot': feet_to_meter = 0.3048 From ca2e23dc42a6307ff6028945b577d008cd3ad563 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 14:07:35 -0700 Subject: [PATCH 10/32] tweaked estimation of resolution. --- tests/test_layer_parameters.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index c939202..49f4bb0 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -109,7 +109,8 @@ def evaluate_resolution__property(obj): data = obj.get_data(BBOX) expected_resolution = doubled_default_resolution - actual_estimated_resolution = estimate_spatial_resolution(data) + actual_estimated_resolution, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) + print (expected_resolution, actual_estimated_resolution, crs_units, diff_distance, y_cells, x_min, y_min, y_max) return expected_resolution, actual_estimated_resolution @@ -153,9 +154,11 @@ def get_class_from_instance(obj): def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) + x_min = None y_min = data['y'].values.min() y_max = data['y'].values.max() + crs_units = None try: crs = CRS.from_string(data.crs) crs_units = crs.axis_info[0].unit_name @@ -164,19 +167,16 @@ def estimate_spatial_resolution(data): crs_units = 'metre' if crs_units == 'metre' or crs_units == 'm': - y_diff = y_max - y_min + diff_distance = y_max - y_min elif crs_units == 'foot': feet_to_meter = 0.3048 - y_diff = (y_max - y_min) * feet_to_meter + diff_distance = (y_max - y_min) / feet_to_meter elif crs_units == 'degree': - lat1 = y_min - lat2 = y_max - lon1 = data['x'].values.min() - lon2 = lon1 - y_diff = get_distance_between_geocoordinates(lat1, lon1, lat2, lon2) + x_min = data['x'].values.min() + diff_distance = get_distance_between_geocoordinates(x_min, y_min, x_min, y_max) else: raise Exception('Unhandled projection units: %s' % crs_units) - ry = round(y_diff / y_cells) + ry = round(diff_distance / y_cells) - return ry + return ry, crs_units, diff_distance, y_cells, x_min, y_min, y_max From 58b2eba1733ea6cf9c720ec318ddeda9ff39d18a Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 15:39:10 -0700 Subject: [PATCH 11/32] improved method for getting crs --- tests/test_layer_parameters.py | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 49f4bb0..3a2809e 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -108,9 +108,11 @@ def evaluate_resolution__property(obj): data = obj.get_data(BBOX) + + expected_resolution = doubled_default_resolution - actual_estimated_resolution, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) - print (expected_resolution, actual_estimated_resolution, crs_units, diff_distance, y_cells, x_min, y_min, y_max) + actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) + print (expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) return expected_resolution, actual_estimated_resolution @@ -158,25 +160,18 @@ def estimate_spatial_resolution(data): y_min = data['y'].values.min() y_max = data['y'].values.max() - crs_units = None - try: - crs = CRS.from_string(data.crs) - crs_units = crs.axis_info[0].unit_name - except: - # if xarray doesn't have crs property, assume meters - crs_units = 'metre' + crs_string = data.rio.crs.data['init'] + crs = CRS.from_string(crs_string) + crs_unit = crs.axis_info[0].unit_name - if crs_units == 'metre' or crs_units == 'm': + if crs_unit == 'metre' or crs_unit == 'm': diff_distance = y_max - y_min - elif crs_units == 'foot': - feet_to_meter = 0.3048 - diff_distance = (y_max - y_min) / feet_to_meter - elif crs_units == 'degree': + elif crs_unit == 'degree': x_min = data['x'].values.min() diff_distance = get_distance_between_geocoordinates(x_min, y_min, x_min, y_max) else: - raise Exception('Unhandled projection units: %s' % crs_units) + raise Exception('Unhandled projection units: %s' % crs_unit) ry = round(diff_distance / y_cells) - return ry, crs_units, diff_distance, y_cells, x_min, y_min, y_max + return ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max From 0da202a5cf5f585393f13fd541c6c168f6254715 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 15:55:05 -0700 Subject: [PATCH 12/32] corrected error in distance calculation --- tests/test_layer_parameters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 3a2809e..f31a89a 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -168,7 +168,7 @@ def estimate_spatial_resolution(data): diff_distance = y_max - y_min elif crs_unit == 'degree': x_min = data['x'].values.min() - diff_distance = get_distance_between_geocoordinates(x_min, y_min, x_min, y_max) + diff_distance = get_distance_between_geocoordinates(y_min, x_min, y_max, x_min) else: raise Exception('Unhandled projection units: %s' % crs_unit) From 134e50bc847ce86b1f61639c534390b418bdf3aa Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 16:09:27 -0700 Subject: [PATCH 13/32] changed determination of crs --- tests/test_layer_parameters.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index f31a89a..402bc24 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -111,8 +111,8 @@ def evaluate_resolution__property(obj): expected_resolution = doubled_default_resolution - actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) - print (expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) + tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) + print (tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) return expected_resolution, actual_estimated_resolution @@ -160,7 +160,10 @@ def estimate_spatial_resolution(data): y_min = data['y'].values.min() y_max = data['y'].values.max() - crs_string = data.rio.crs.data['init'] + try: + crs_string = data.crs + except: + crs_string = data.rio.crs.data['init'] crs = CRS.from_string(crs_string) crs_unit = crs.axis_info[0].unit_name @@ -174,4 +177,5 @@ def estimate_spatial_resolution(data): ry = round(diff_distance / y_cells) - return ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max + tt = type(data) + return tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max From 8f1cbbe9a3c7e12de919c42da4519dd1ea667d63 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 17:42:04 -0700 Subject: [PATCH 14/32] reduced to natural areas --- tests/test_layer_parameters.py | 126 +++++++++++++++++---------------- 1 file changed, 64 insertions(+), 62 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 402bc24..abcc510 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -30,70 +30,70 @@ COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 -def test_albedo_spatial_resolution(): - class_instance = Albedo() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_alos_dsm_spatial_resolution(): - class_instance = AlosDSM() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_average_net_building_height_spatial_resolution(): - class_instance = AverageNetBuildingHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_built_up_height_spatial_resolution(): - class_instance = BuiltUpHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_esa_world_cover_spatial_resolution(): - class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_land_surface_temperature_spatial_resolution(): - class_instance = LandSurfaceTemperature() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_nasa_dem_spatial_resolution(): - class_instance = NasaDEM() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - +# def test_albedo_spatial_resolution(): +# class_instance = Albedo() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_alos_dsm_spatial_resolution(): +# class_instance = AlosDSM() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_average_net_building_height_spatial_resolution(): +# class_instance = AverageNetBuildingHeight() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_built_up_height_spatial_resolution(): +# class_instance = BuiltUpHeight() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_esa_world_cover_spatial_resolution(): +# class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_land_surface_temperature_spatial_resolution(): +# class_instance = LandSurfaceTemperature() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_nasa_dem_spatial_resolution(): +# class_instance = NasaDEM() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# def test_natural_areas_spatial_resolution(): class_instance = NaturalAreas() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert doubled_default_resolution == actual_estimated_resolution -def test_ndvi_sentinel2_spatial_resolution(): - class_instance = NdviSentinel2(year=2023) - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_tree_canopy_height_spatial_resolution(): - class_instance = TreeCanopyHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_tree_cover_spatial_resolution(): - class_instance = TreeCover() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_urban_land_use_spatial_resolution(): - class_instance = UrbanLandUse() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution - -def test_world_pop_spatial_resolution(): - class_instance = WorldPop() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution +# def test_ndvi_sentinel2_spatial_resolution(): +# class_instance = NdviSentinel2(year=2023) +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_tree_canopy_height_spatial_resolution(): +# class_instance = TreeCanopyHeight() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_tree_cover_spatial_resolution(): +# class_instance = TreeCover() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_urban_land_use_spatial_resolution(): +# class_instance = UrbanLandUse() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution +# +# def test_world_pop_spatial_resolution(): +# class_instance = WorldPop() +# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) +# assert doubled_default_resolution == actual_estimated_resolution def evaluate_resolution__property(obj): is_valid, except_str = validate_layer_instance(obj) @@ -111,8 +111,8 @@ def evaluate_resolution__property(obj): expected_resolution = doubled_default_resolution - tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) - print (tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) + method, tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) + print (method, tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) return expected_resolution, actual_estimated_resolution @@ -161,8 +161,10 @@ def estimate_spatial_resolution(data): y_max = data['y'].values.max() try: + method = 'a' crs_string = data.crs except: + method = 'b' crs_string = data.rio.crs.data['init'] crs = CRS.from_string(crs_string) crs_unit = crs.axis_info[0].unit_name @@ -178,4 +180,4 @@ def estimate_spatial_resolution(data): ry = round(diff_distance / y_cells) tt = type(data) - return tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max + return method, tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max From 4fdc35cc8bf8a23db93af5c7d57f91cd99925fae Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 18:11:03 -0700 Subject: [PATCH 15/32] enabled all tests --- tests/test_layer_parameters.py | 118 ++++++++++++++++----------------- 1 file changed, 59 insertions(+), 59 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index abcc510..12d4ae5 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -30,70 +30,70 @@ COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 -# def test_albedo_spatial_resolution(): -# class_instance = Albedo() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_alos_dsm_spatial_resolution(): -# class_instance = AlosDSM() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_average_net_building_height_spatial_resolution(): -# class_instance = AverageNetBuildingHeight() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_built_up_height_spatial_resolution(): -# class_instance = BuiltUpHeight() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_esa_world_cover_spatial_resolution(): -# class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_land_surface_temperature_spatial_resolution(): -# class_instance = LandSurfaceTemperature() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_nasa_dem_spatial_resolution(): -# class_instance = NasaDEM() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# +def test_albedo_spatial_resolution(): + class_instance = Albedo() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_alos_dsm_spatial_resolution(): + class_instance = AlosDSM() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_average_net_building_height_spatial_resolution(): + class_instance = AverageNetBuildingHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_built_up_height_spatial_resolution(): + class_instance = BuiltUpHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_esa_world_cover_spatial_resolution(): + class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_land_surface_temperature_spatial_resolution(): + class_instance = LandSurfaceTemperature() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_nasa_dem_spatial_resolution(): + class_instance = NasaDEM() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + def test_natural_areas_spatial_resolution(): class_instance = NaturalAreas() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert doubled_default_resolution == actual_estimated_resolution -# def test_ndvi_sentinel2_spatial_resolution(): -# class_instance = NdviSentinel2(year=2023) -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_tree_canopy_height_spatial_resolution(): -# class_instance = TreeCanopyHeight() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_tree_cover_spatial_resolution(): -# class_instance = TreeCover() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_urban_land_use_spatial_resolution(): -# class_instance = UrbanLandUse() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution -# -# def test_world_pop_spatial_resolution(): -# class_instance = WorldPop() -# doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) -# assert doubled_default_resolution == actual_estimated_resolution +def test_ndvi_sentinel2_spatial_resolution(): + class_instance = NdviSentinel2(year=2023) + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_tree_canopy_height_spatial_resolution(): + class_instance = TreeCanopyHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_tree_cover_spatial_resolution(): + class_instance = TreeCover() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_urban_land_use_spatial_resolution(): + class_instance = UrbanLandUse() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution + +def test_world_pop_spatial_resolution(): + class_instance = WorldPop() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert doubled_default_resolution == actual_estimated_resolution def evaluate_resolution__property(obj): is_valid, except_str = validate_layer_instance(obj) From 9f569c9ced51d34c7763e3373255a04bb572fc65 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 21:30:53 -0700 Subject: [PATCH 16/32] changed determination of crs_units --- tests/resources/bbox_constants.py | 12 ++++++++---- tests/test_layer_parameters.py | 29 +++++++++++++++++------------ 2 files changed, 25 insertions(+), 16 deletions(-) diff --git a/tests/resources/bbox_constants.py b/tests/resources/bbox_constants.py index 9aaf8ad..789ab48 100644 --- a/tests/resources/bbox_constants.py +++ b/tests/resources/bbox_constants.py @@ -9,10 +9,14 @@ ) BBOX_BRA_SALVADOR_ADM4 = ( - -38.647320153390055, - -13.01748678217598787, - -38.3041637148564007, - -12.75607703449720631 + -38.647320153390055,-13.01748678217598787, + -38.3041637148564007,-12.75607703449720631 +) + +# UTM Zones 22S and 23S +BBOX_BRA_BRASILIA = ( + -48.07651,-15.89788 + -47.83736,-15.71919 ) BBOX_SMALL_TEST = ( diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 12d4ae5..b831862 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -15,7 +15,7 @@ UrbanLandUse, WorldPop, OpenStreetMap ) -from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 +from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1, BBOX_BRA_SALVADOR_ADM4, BBOX_BRA_BRASILIA from tests.tools.spatial_tools import get_distance_between_geocoordinates """ @@ -29,6 +29,7 @@ COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 +# BBOX = BBOX_BRA_BRASILIA def test_albedo_spatial_resolution(): class_instance = Albedo() @@ -111,8 +112,8 @@ def evaluate_resolution__property(obj): expected_resolution = doubled_default_resolution - method, tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) - print (method, tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) + tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) + print (tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) return expected_resolution, actual_estimated_resolution @@ -160,14 +161,18 @@ def estimate_spatial_resolution(data): y_min = data['y'].values.min() y_max = data['y'].values.max() - try: - method = 'a' - crs_string = data.crs - except: - method = 'b' - crs_string = data.rio.crs.data['init'] - crs = CRS.from_string(crs_string) - crs_unit = crs.axis_info[0].unit_name + crs = data.rio.crs + crs_unit = data.rio.crs.linear_units + # try: + # method = 'a' + # # crs_string = data.crs + # crs = data.rio.crs + # crs_unit = data.rio.crs.linear_units + # except: + # method = 'b' + # crs_string = data.rio.crs.data['init'] + # crs = CRS.from_string(crs_string) + # crs_unit = crs.axis_info[0].unit_name if crs_unit == 'metre' or crs_unit == 'm': diff_distance = y_max - y_min @@ -180,4 +185,4 @@ def estimate_spatial_resolution(data): ry = round(diff_distance / y_cells) tt = type(data) - return method, tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max + return tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max From 2168ef8ac72f059c4d5d1035ef46eb98bdca2496 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 22:13:46 -0700 Subject: [PATCH 17/32] changed computation of coordinates --- tests/test_layer_parameters.py | 36 +++++++++++++++++----------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index b831862..fad870f 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -112,8 +112,8 @@ def evaluate_resolution__property(obj): expected_resolution = doubled_default_resolution - tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) - print (tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) + method, tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) + print (method, tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) return expected_resolution, actual_estimated_resolution @@ -158,21 +158,21 @@ def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) x_min = None - y_min = data['y'].values.min() - y_max = data['y'].values.max() - - crs = data.rio.crs - crs_unit = data.rio.crs.linear_units - # try: - # method = 'a' - # # crs_string = data.crs - # crs = data.rio.crs - # crs_unit = data.rio.crs.linear_units - # except: - # method = 'b' - # crs_string = data.rio.crs.data['init'] - # crs = CRS.from_string(crs_string) - # crs_unit = crs.axis_info[0].unit_name + # y_min = data['y'].values.min() + # y_max = data['y'].values.max() + y_min = data.coords['y'].values.min() + y_max = data.coords['y'].values.max() + + try: + method = 'a' + crs_string = data.crs + crs = CRS.from_string(crs_string) + crs_unit = crs.axis_info[0].unit_name + except: + method = 'b' + # crs_string = data.rio.crs.data['init'] + crs = data.rio.crs + crs_unit = data.rio.crs.linear_units if crs_unit == 'metre' or crs_unit == 'm': diff_distance = y_max - y_min @@ -185,4 +185,4 @@ def estimate_spatial_resolution(data): ry = round(diff_distance / y_cells) tt = type(data) - return tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max + return method, tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max From 336a5a5967997a326619d6a21c7dc7c385b6a0d9 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 22:28:52 -0700 Subject: [PATCH 18/32] changed handling of 'm' unit --- tests/test_layer_parameters.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index fad870f..949b192 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,3 +1,4 @@ +import pytest from pyproj import CRS from city_metrix.layers import ( Layer, @@ -29,72 +30,72 @@ COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 -# BBOX = BBOX_BRA_BRASILIA +RESOLUTION_TOLERANCE = 1 def test_albedo_spatial_resolution(): class_instance = Albedo() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_alos_dsm_spatial_resolution(): class_instance = AlosDSM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_average_net_building_height_spatial_resolution(): class_instance = AverageNetBuildingHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_built_up_height_spatial_resolution(): class_instance = BuiltUpHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_esa_world_cover_spatial_resolution(): class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_land_surface_temperature_spatial_resolution(): class_instance = LandSurfaceTemperature() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_nasa_dem_spatial_resolution(): class_instance = NasaDEM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_natural_areas_spatial_resolution(): class_instance = NaturalAreas() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_ndvi_sentinel2_spatial_resolution(): class_instance = NdviSentinel2(year=2023) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_tree_canopy_height_spatial_resolution(): class_instance = TreeCanopyHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_tree_cover_spatial_resolution(): class_instance = TreeCover() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_urban_land_use_spatial_resolution(): class_instance = UrbanLandUse() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_world_pop_spatial_resolution(): class_instance = WorldPop() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) - assert doubled_default_resolution == actual_estimated_resolution + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def evaluate_resolution__property(obj): is_valid, except_str = validate_layer_instance(obj) @@ -174,9 +175,9 @@ def estimate_spatial_resolution(data): crs = data.rio.crs crs_unit = data.rio.crs.linear_units - if crs_unit == 'metre' or crs_unit == 'm': + if crs_unit == 'metre': diff_distance = y_max - y_min - elif crs_unit == 'degree': + elif crs_unit == 'degree' or crs_unit == 'm': x_min = data['x'].values.min() diff_distance = get_distance_between_geocoordinates(y_min, x_min, y_max, x_min) else: From 976055759ef395563146f88002f5bde80fad1352 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 22:45:42 -0700 Subject: [PATCH 19/32] changed determination of crs_unit --- tests/test_layer_parameters.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 949b192..a0356a6 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -171,14 +171,17 @@ def estimate_spatial_resolution(data): crs_unit = crs.axis_info[0].unit_name except: method = 'b' - # crs_string = data.rio.crs.data['init'] - crs = data.rio.crs - crs_unit = data.rio.crs.linear_units + crs_string = data.rio.crs.data['init'] + # crs = data.rio.crs + # crs_unit = data.rio.crs.linear_units + crs = CRS.from_string(crs_string) + crs_unit = crs.axis_info[0].unit_name if crs_unit == 'metre': diff_distance = y_max - y_min elif crs_unit == 'degree' or crs_unit == 'm': - x_min = data['x'].values.min() + # x_min = data['x'].values.min() + x_min = data.coords['x'].values.min() diff_distance = get_distance_between_geocoordinates(y_min, x_min, y_max, x_min) else: raise Exception('Unhandled projection units: %s' % crs_unit) From ac80767b6a0cdf64b3b993e9f0b7addcb2c3ffd6 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 22:57:57 -0700 Subject: [PATCH 20/32] forced to fail --- tests/test_layer_parameters.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index a0356a6..4615747 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -36,66 +36,79 @@ def test_albedo_spatial_resolution(): class_instance = Albedo() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_alos_dsm_spatial_resolution(): class_instance = AlosDSM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_average_net_building_height_spatial_resolution(): class_instance = AverageNetBuildingHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_built_up_height_spatial_resolution(): class_instance = BuiltUpHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_esa_world_cover_spatial_resolution(): class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_land_surface_temperature_spatial_resolution(): class_instance = LandSurfaceTemperature() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_nasa_dem_spatial_resolution(): class_instance = NasaDEM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_natural_areas_spatial_resolution(): class_instance = NaturalAreas() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_ndvi_sentinel2_spatial_resolution(): class_instance = NdviSentinel2(year=2023) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_tree_canopy_height_spatial_resolution(): class_instance = TreeCanopyHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_tree_cover_spatial_resolution(): class_instance = TreeCover() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_urban_land_use_spatial_resolution(): class_instance = UrbanLandUse() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def test_world_pop_spatial_resolution(): class_instance = WorldPop() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + assert 1==2 def evaluate_resolution__property(obj): is_valid, except_str = validate_layer_instance(obj) @@ -159,8 +172,6 @@ def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) x_min = None - # y_min = data['y'].values.min() - # y_max = data['y'].values.max() y_min = data.coords['y'].values.min() y_max = data.coords['y'].values.max() @@ -172,15 +183,12 @@ def estimate_spatial_resolution(data): except: method = 'b' crs_string = data.rio.crs.data['init'] - # crs = data.rio.crs - # crs_unit = data.rio.crs.linear_units crs = CRS.from_string(crs_string) crs_unit = crs.axis_info[0].unit_name if crs_unit == 'metre': diff_distance = y_max - y_min elif crs_unit == 'degree' or crs_unit == 'm': - # x_min = data['x'].values.min() x_min = data.coords['x'].values.min() diff_distance = get_distance_between_geocoordinates(y_min, x_min, y_max, x_min) else: From e7956cdc52867657095947daddd708da9206b990 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Fri, 30 Aug 2024 23:18:17 -0700 Subject: [PATCH 21/32] try to identify wkt projection --- tests/test_layer_parameters.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 4615747..a4b7763 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -178,7 +178,10 @@ def estimate_spatial_resolution(data): try: method = 'a' crs_string = data.crs - crs = CRS.from_string(crs_string) + if crs_string.startswith('PROJCS['): + crs = CRS.from_wkt(crs_string) + else: + crs = CRS.from_string(crs_string) crs_unit = crs.axis_info[0].unit_name except: method = 'b' From 0aefba2f0beb717d948dc17791f4903841c26f68 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Tue, 3 Sep 2024 10:22:33 -0700 Subject: [PATCH 22/32] trying to isolate projection issue --- city_metrix/layers/layer.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 01ad6e4..74eedfe 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -275,8 +275,18 @@ def get_image_collection( :return: """ +# TODO may not be reprojecting!!!! + crs = get_utm_zone_epsg(bbox) + ds_orig = xr.open_dataset( + image_collection, + engine='ee', + scale=scale, + geometry=ee.Geometry.Rectangle(*bbox), + chunks={'X': 512, 'Y': 512}, + ) + ds = xr.open_dataset( image_collection, engine='ee', @@ -290,6 +300,15 @@ def get_image_collection( print(f"Extracting layer {name} from Google Earth Engine for bbox {bbox}:") data = ds.compute() + x_val_orig = ds_orig.coords['lon'].values.min() + x_val_modified = ds.coords['X'].values.min() + + crs_orig = ds_orig.attrs['crs'] + crs_modified = ds.attrs['crs'] + print('CRS: %s, Org_x: %s' % (crs_orig, x_val_orig)) + print('CRS: %s, Modified_x: %s' % (crs_modified, x_val_modified)) + + # get in rioxarray format data = data.squeeze("time").transpose("Y", "X").rename({'X': 'x', 'Y': 'y'}) From 34759aa0c0ab9b65920544417a0251a158434ec5 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Tue, 3 Sep 2024 10:33:23 -0700 Subject: [PATCH 23/32] more attempts to isolate projection issue --- city_metrix/layers/layer.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 74eedfe..7678fb6 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -300,13 +300,16 @@ def get_image_collection( print(f"Extracting layer {name} from Google Earth Engine for bbox {bbox}:") data = ds.compute() + crs_orig = ds_orig.attrs['crs'] + crs_modified = ds.attrs['crs'] + print('CRS: %s' % (crs_orig)) + print('CRS: %s' % (crs_modified)) + x_val_orig = ds_orig.coords['lon'].values.min() x_val_modified = ds.coords['X'].values.min() + print('Org_x: %s' % (x_val_orig)) + print('Modified_x: %s' % (x_val_modified)) - crs_orig = ds_orig.attrs['crs'] - crs_modified = ds.attrs['crs'] - print('CRS: %s, Org_x: %s' % (crs_orig, x_val_orig)) - print('CRS: %s, Modified_x: %s' % (crs_modified, x_val_modified)) # get in rioxarray format From 6aab4c2856421e545016a31418d68b467d9fdd37 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Tue, 3 Sep 2024 11:27:23 -0700 Subject: [PATCH 24/32] another attempts to isolate projection issue --- city_metrix/layers/layer.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 7678fb6..2e5a53e 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -284,7 +284,6 @@ def get_image_collection( engine='ee', scale=scale, geometry=ee.Geometry.Rectangle(*bbox), - chunks={'X': 512, 'Y': 512}, ) ds = xr.open_dataset( @@ -305,10 +304,10 @@ def get_image_collection( print('CRS: %s' % (crs_orig)) print('CRS: %s' % (crs_modified)) - x_val_orig = ds_orig.coords['lon'].values.min() - x_val_modified = ds.coords['X'].values.min() - print('Org_x: %s' % (x_val_orig)) - print('Modified_x: %s' % (x_val_modified)) + # x_val_orig = ds_orig.coords['lon'].values.min() + # x_val_modified = ds.coords['X'].values.min() + # print('Org_x: %s' % (x_val_orig)) + # print('Modified_x: %s' % (x_val_modified)) From 7b2bbe6225ee3b8106ae4faf3b683bc8fb20f327 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Tue, 3 Sep 2024 21:28:21 -0700 Subject: [PATCH 25/32] Identified projection error as xee version --- .github/requirements.txt | 4 +-- city_metrix/layers/layer.py | 24 ++-------------- tests/test_layer_parameters.py | 51 +++++++++++++++------------------- 3 files changed, 26 insertions(+), 53 deletions(-) diff --git a/.github/requirements.txt b/.github/requirements.txt index 9928a46..171974e 100644 --- a/.github/requirements.txt +++ b/.github/requirements.txt @@ -1,12 +1,12 @@ earthengine-api==0.1.408 geocube==0.4.2 -geopandas==0.14.1 +geopandas==0.14.4 rioxarray==0.15.0 odc-stac==0.3.8 pystac-client==0.7.5 pytest==7.4.3 xarray-spatial==0.3.7 -xee==0.0.3 +xee==0.0.15 utm==0.7.0 osmnx==1.9.3 dask[complete]==2023.11.0 diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 2e5a53e..247fb87 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -17,6 +17,7 @@ import utm import shapely.geometry as geometry import pandas as pd +from pyproj import Transformer MAX_TILE_SIZE = 0.5 @@ -274,18 +275,9 @@ def get_image_collection( :param name: optional name to print while reporting progress :return: """ - -# TODO may not be reprojecting!!!! - crs = get_utm_zone_epsg(bbox) - ds_orig = xr.open_dataset( - image_collection, - engine='ee', - scale=scale, - geometry=ee.Geometry.Rectangle(*bbox), - ) - + # Use a higher version for xee https://github.com/google/Xee/issues/118 ds = xr.open_dataset( image_collection, engine='ee', @@ -299,18 +291,6 @@ def get_image_collection( print(f"Extracting layer {name} from Google Earth Engine for bbox {bbox}:") data = ds.compute() - crs_orig = ds_orig.attrs['crs'] - crs_modified = ds.attrs['crs'] - print('CRS: %s' % (crs_orig)) - print('CRS: %s' % (crs_modified)) - - # x_val_orig = ds_orig.coords['lon'].values.min() - # x_val_modified = ds.coords['X'].values.min() - # print('Org_x: %s' % (x_val_orig)) - # print('Modified_x: %s' % (x_val_modified)) - - - # get in rioxarray format data = data.squeeze("time").transpose("Y", "X").rename({'X': 'x', 'Y': 'y'}) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index a4b7763..2710a0d 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -36,98 +36,98 @@ def test_albedo_spatial_resolution(): class_instance = Albedo() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_alos_dsm_spatial_resolution(): class_instance = AlosDSM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_average_net_building_height_spatial_resolution(): class_instance = AverageNetBuildingHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_built_up_height_spatial_resolution(): class_instance = BuiltUpHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_esa_world_cover_spatial_resolution(): class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_land_surface_temperature_spatial_resolution(): class_instance = LandSurfaceTemperature() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_nasa_dem_spatial_resolution(): class_instance = NasaDEM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_natural_areas_spatial_resolution(): class_instance = NaturalAreas() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_ndvi_sentinel2_spatial_resolution(): class_instance = NdviSentinel2(year=2023) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_tree_canopy_height_spatial_resolution(): class_instance = TreeCanopyHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_tree_cover_spatial_resolution(): class_instance = TreeCover() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_urban_land_use_spatial_resolution(): class_instance = UrbanLandUse() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 def test_world_pop_spatial_resolution(): class_instance = WorldPop() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - assert 1==2 + # assert 1==2 + +# def evaluate_resolution__property(obj): +# data = obj.get_data(BBOX) def evaluate_resolution__property(obj): is_valid, except_str = validate_layer_instance(obj) if is_valid is False: raise Exception(except_str) + # Double the default scale for testing cls = get_class_from_instance(obj) - - # Double the spatial_resolution for the specified Class doubled_default_resolution = 2 * cls.spatial_resolution obj.spatial_resolution=doubled_default_resolution data = obj.get_data(BBOX) - - expected_resolution = doubled_default_resolution - method, tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max = estimate_spatial_resolution(data) - print (method, tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, x_min, y_min, y_max) + tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, y_min, y_max = estimate_spatial_resolution(data) + print (tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, y_min, y_max) return expected_resolution, actual_estimated_resolution @@ -170,34 +170,27 @@ def get_class_from_instance(obj): def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) - - x_min = None y_min = data.coords['y'].values.min() y_max = data.coords['y'].values.max() try: - method = 'a' crs_string = data.crs if crs_string.startswith('PROJCS['): crs = CRS.from_wkt(crs_string) else: crs = CRS.from_string(crs_string) - crs_unit = crs.axis_info[0].unit_name + crs_unit = crs.axis_info[0].unit_name except: - method = 'b' crs_string = data.rio.crs.data['init'] crs = CRS.from_string(crs_string) crs_unit = crs.axis_info[0].unit_name if crs_unit == 'metre': diff_distance = y_max - y_min - elif crs_unit == 'degree' or crs_unit == 'm': - x_min = data.coords['x'].values.min() - diff_distance = get_distance_between_geocoordinates(y_min, x_min, y_max, x_min) else: - raise Exception('Unhandled projection units: %s' % crs_unit) + raise Exception('Unhandled projection units: %s for projection: %s' % (crs_unit, crs_string)) ry = round(diff_distance / y_cells) tt = type(data) - return method, tt, ry, crs, crs_unit, diff_distance, y_cells, x_min, y_min, y_max + return tt, ry, crs, crs_unit, diff_distance, y_cells, y_min, y_max From 409e706c73395e9717f2b131e2281a75429fa39e Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Tue, 3 Sep 2024 22:14:26 -0700 Subject: [PATCH 26/32] Cleaned up code and updated requirements --- city_metrix/layers/layer.py | 4 ++-- environment.yml | 4 ++-- tests/test_layer_parameters.py | 43 +++++++--------------------------- tests/tools/spatial_tools.py | 16 ------------- 4 files changed, 12 insertions(+), 55 deletions(-) delete mode 100644 tests/tools/spatial_tools.py diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 247fb87..86590c8 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -17,7 +17,6 @@ import utm import shapely.geometry as geometry import pandas as pd -from pyproj import Transformer MAX_TILE_SIZE = 0.5 @@ -275,9 +274,10 @@ def get_image_collection( :param name: optional name to print while reporting progress :return: """ + crs = get_utm_zone_epsg(bbox) - # Use a higher version for xee https://github.com/google/Xee/issues/118 + # See link regarding bug in crs specification shttps://github.com/google/Xee/issues/118 ds = xr.open_dataset( image_collection, engine='ee', diff --git a/environment.yml b/environment.yml index 24ec040..e68a8fc 100644 --- a/environment.yml +++ b/environment.yml @@ -5,13 +5,13 @@ dependencies: - python=3.10 - earthengine-api=0.1.379 - geocube=0.4.2 - - geopandas=0.14.1 + - geopandas=0.14.4 - rioxarray=0.15.0 - odc-stac=0.3.8 - pystac-client=0.7.5 - pytest=7.4.3 - xarray-spatial=0.3.7 - - xee=0.0.3 + - xee=0.0.15 - utm=0.7.0 - osmnx=1.9.0 - dask[complete]=2023.11.0 diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 2710a0d..2238bad 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -16,8 +16,7 @@ UrbanLandUse, WorldPop, OpenStreetMap ) -from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1, BBOX_BRA_SALVADOR_ADM4, BBOX_BRA_BRASILIA -from tests.tools.spatial_tools import get_distance_between_geocoordinates +from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 """ Evaluation of spatial_resolution property @@ -36,82 +35,66 @@ def test_albedo_spatial_resolution(): class_instance = Albedo() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_alos_dsm_spatial_resolution(): class_instance = AlosDSM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_average_net_building_height_spatial_resolution(): class_instance = AverageNetBuildingHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_built_up_height_spatial_resolution(): class_instance = BuiltUpHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_esa_world_cover_spatial_resolution(): class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_land_surface_temperature_spatial_resolution(): class_instance = LandSurfaceTemperature() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_nasa_dem_spatial_resolution(): class_instance = NasaDEM() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_natural_areas_spatial_resolution(): class_instance = NaturalAreas() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_ndvi_sentinel2_spatial_resolution(): class_instance = NdviSentinel2(year=2023) doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_tree_canopy_height_spatial_resolution(): class_instance = TreeCanopyHeight() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_tree_cover_spatial_resolution(): class_instance = TreeCover() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_urban_land_use_spatial_resolution(): class_instance = UrbanLandUse() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 def test_world_pop_spatial_resolution(): class_instance = WorldPop() doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution - # assert 1==2 - -# def evaluate_resolution__property(obj): -# data = obj.get_data(BBOX) def evaluate_resolution__property(obj): is_valid, except_str = validate_layer_instance(obj) @@ -126,10 +109,9 @@ def evaluate_resolution__property(obj): data = obj.get_data(BBOX) expected_resolution = doubled_default_resolution - tt, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, y_min, y_max = estimate_spatial_resolution(data) - print (tt, expected_resolution, actual_estimated_resolution, crs, crs_units, diff_distance, y_cells, y_min, y_max) + estimated_actual_resolution = estimate_spatial_resolution(data) - return expected_resolution, actual_estimated_resolution + return expected_resolution, estimated_actual_resolution def test_function_validate_layer_instance(): @@ -173,24 +155,15 @@ def estimate_spatial_resolution(data): y_min = data.coords['y'].values.min() y_max = data.coords['y'].values.max() - try: - crs_string = data.crs - if crs_string.startswith('PROJCS['): - crs = CRS.from_wkt(crs_string) - else: - crs = CRS.from_string(crs_string) - crs_unit = crs.axis_info[0].unit_name - except: - crs_string = data.rio.crs.data['init'] - crs = CRS.from_string(crs_string) - crs_unit = crs.axis_info[0].unit_name + crs_string = data.rio.crs.data['init'] + crs = CRS.from_string(crs_string) + crs_unit = crs.axis_info[0].unit_name if crs_unit == 'metre': diff_distance = y_max - y_min else: raise Exception('Unhandled projection units: %s for projection: %s' % (crs_unit, crs_string)) - ry = round(diff_distance / y_cells) + estimated_actual_resolution = round(diff_distance / y_cells) - tt = type(data) - return tt, ry, crs, crs_unit, diff_distance, y_cells, y_min, y_max + return estimated_actual_resolution diff --git a/tests/tools/spatial_tools.py b/tests/tools/spatial_tools.py deleted file mode 100644 index 73ab093..0000000 --- a/tests/tools/spatial_tools.py +++ /dev/null @@ -1,16 +0,0 @@ -import math - -EARTH_RADIUS = 6378.137 # Radius of earth in KM - -def get_distance_between_geocoordinates(lat1, lon1, lat2, lon2): - dLat = lat2 * math.pi / 180 - lat1 * math.pi / 180 - dLon = lon2 * math.pi / 180 - lon1 * math.pi / 180 - a = math.sin(dLat/2) * math.sin(dLat/2) + \ - math.cos(lat1 * math.pi / 180) * \ - math.cos(lat2 * math.pi / 180) * \ - math.sin(dLon/2) * math.sin(dLon/2) - c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) - d = EARTH_RADIUS * c - return d * 1000 # meters - - From 85b02f401753c6468ad5a2661b4fe7225eb6e6ad Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Tue, 3 Sep 2024 22:32:08 -0700 Subject: [PATCH 27/32] Squashing commits related to understanding GH Actions --- .github/requirements.txt | 6 +- .github/workflows/dev_ci_cd.yml | 4 +- city_metrix/layers/layer.py | 1 + environment.yml | 6 +- setup.py | 2 +- tests/resources/bbox_constants.py | 12 +- tests/test_layer_parameters.py | 215 +++++++++++++++++------------- tests/tools/spatial_tools.py | 16 --- 8 files changed, 142 insertions(+), 120 deletions(-) delete mode 100644 tests/tools/spatial_tools.py diff --git a/.github/requirements.txt b/.github/requirements.txt index 9928a46..15d8543 100644 --- a/.github/requirements.txt +++ b/.github/requirements.txt @@ -1,12 +1,12 @@ earthengine-api==0.1.408 geocube==0.4.2 -geopandas==0.14.1 +geopandas==0.14.4 rioxarray==0.15.0 odc-stac==0.3.8 pystac-client==0.7.5 pytest==7.4.3 xarray-spatial==0.3.7 -xee==0.0.3 +xee==0.0.15 utm==0.7.0 osmnx==1.9.3 dask[complete]==2023.11.0 @@ -17,4 +17,4 @@ pip==23.3.1 boto3==1.34.124 scikit-learn==1.5.0 overturemaps==0.6.0 -git+https://github.com/isciences/exactextract \ No newline at end of file +exactextract==0.2.0.dev252 diff --git a/.github/workflows/dev_ci_cd.yml b/.github/workflows/dev_ci_cd.yml index 6d55444..c57fd39 100644 --- a/.github/workflows/dev_ci_cd.yml +++ b/.github/workflows/dev_ci_cd.yml @@ -14,9 +14,9 @@ jobs: python-version: ["3.10"] steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v3 + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} - name: Install Linux dependencies diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 01ad6e4..86590c8 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -277,6 +277,7 @@ def get_image_collection( crs = get_utm_zone_epsg(bbox) + # See link regarding bug in crs specification shttps://github.com/google/Xee/issues/118 ds = xr.open_dataset( image_collection, engine='ee', diff --git a/environment.yml b/environment.yml index 24ec040..2434184 100644 --- a/environment.yml +++ b/environment.yml @@ -5,13 +5,13 @@ dependencies: - python=3.10 - earthengine-api=0.1.379 - geocube=0.4.2 - - geopandas=0.14.1 + - geopandas=0.14.4 - rioxarray=0.15.0 - odc-stac=0.3.8 - pystac-client=0.7.5 - pytest=7.4.3 - xarray-spatial=0.3.7 - - xee=0.0.3 + - xee=0.0.15 - utm=0.7.0 - osmnx=1.9.0 - dask[complete]=2023.11.0 @@ -22,6 +22,6 @@ dependencies: - pip=23.3.1 - boto3=1.34.124 - scikit-learn=1.5.0 + - exactextract=0.2.0.dev252 - pip: - - git+https://github.com/isciences/exactextract - overturemaps==0.6.0 diff --git a/setup.py b/setup.py index 3124621..179b595 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ "s3fs", "dask>=2023.11.0", "boto3", - "exactextract", + "exactextract<=0.2.0.dev252", "overturemaps", "scikit-learn>=1.5.0", ], diff --git a/tests/resources/bbox_constants.py b/tests/resources/bbox_constants.py index 9aaf8ad..789ab48 100644 --- a/tests/resources/bbox_constants.py +++ b/tests/resources/bbox_constants.py @@ -9,10 +9,14 @@ ) BBOX_BRA_SALVADOR_ADM4 = ( - -38.647320153390055, - -13.01748678217598787, - -38.3041637148564007, - -12.75607703449720631 + -38.647320153390055,-13.01748678217598787, + -38.3041637148564007,-12.75607703449720631 +) + +# UTM Zones 22S and 23S +BBOX_BRA_BRASILIA = ( + -48.07651,-15.89788 + -47.83736,-15.71919 ) BBOX_SMALL_TEST = ( diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 043d917..2238bad 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,3 +1,4 @@ +import pytest from pyproj import CRS from city_metrix.layers import ( Layer, @@ -5,8 +6,7 @@ AlosDSM, AverageNetBuildingHeight, BuiltUpHeight, - EsaWorldCover, - EsaWorldCoverClass, + EsaWorldCover, EsaWorldCoverClass, LandSurfaceTemperature, NasaDEM, NaturalAreas, @@ -14,123 +14,156 @@ TreeCanopyHeight, TreeCover, UrbanLandUse, - WorldPop + WorldPop, OpenStreetMap ) from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 -from tests.tools.spatial_tools import get_distance_between_geocoordinates """ -Note: To add a test for another scalable layer that has the spatial_resolution property: +Evaluation of spatial_resolution property +To add a test for a scalable layer that has the spatial_resolution property: 1. Add the class name to the city_metrix.layers import statement above -2. Specify a minimal class instance in the set below. Do no specify the spatial_resolution - property in the instance definition. +2. Copy an existing test_*_spatial_resolution() test + a. rename for the new layer + b. specify a minimal class instance for the layer, not specifying the spatial_resolution attribute. """ -CLASSES_WITH_spatial_resolution_PROPERTY = \ - { - # 'Albedo()', - # 'AlosDSM()', - # 'AverageNetBuildingHeight()', - # 'BuiltUpHeight()', - 'EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)', - # 'LandSurfaceTemperature()', - # 'NasaDEM()', - # 'NaturalAreas()', - # 'NdviSentinel2(year=2023)', - # 'TreeCanopyHeight()', - # 'TreeCover()', - # 'UrbanLandUse()', - # 'WorldPop()' - } COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 - -def test_spatial_resolution_for_all_scalable_layers(): - for class_instance_str in CLASSES_WITH_spatial_resolution_PROPERTY: - is_valid, except_str = validate_layer_instance(class_instance_str) - if is_valid is False: - raise Exception(except_str) - - class_instance = eval(class_instance_str) - - # Double the spatial_resolution for the specified Class - doubled_default_resolution = 2 * class_instance.spatial_resolution - class_instance.spatial_resolution=doubled_default_resolution - - evaluate_layer(class_instance, doubled_default_resolution) +RESOLUTION_TOLERANCE = 1 + +def test_albedo_spatial_resolution(): + class_instance = Albedo() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_alos_dsm_spatial_resolution(): + class_instance = AlosDSM() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_average_net_building_height_spatial_resolution(): + class_instance = AverageNetBuildingHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_built_up_height_spatial_resolution(): + class_instance = BuiltUpHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_esa_world_cover_spatial_resolution(): + class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_land_surface_temperature_spatial_resolution(): + class_instance = LandSurfaceTemperature() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_nasa_dem_spatial_resolution(): + class_instance = NasaDEM() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_natural_areas_spatial_resolution(): + class_instance = NaturalAreas() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_ndvi_sentinel2_spatial_resolution(): + class_instance = NdviSentinel2(year=2023) + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_tree_canopy_height_spatial_resolution(): + class_instance = TreeCanopyHeight() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_tree_cover_spatial_resolution(): + class_instance = TreeCover() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_urban_land_use_spatial_resolution(): + class_instance = UrbanLandUse() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_world_pop_spatial_resolution(): + class_instance = WorldPop() + doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def evaluate_resolution__property(obj): + is_valid, except_str = validate_layer_instance(obj) + if is_valid is False: + raise Exception(except_str) + + # Double the default scale for testing + cls = get_class_from_instance(obj) + doubled_default_resolution = 2 * cls.spatial_resolution + obj.spatial_resolution=doubled_default_resolution + + data = obj.get_data(BBOX) + + expected_resolution = doubled_default_resolution + estimated_actual_resolution = estimate_spatial_resolution(data) + + return expected_resolution, estimated_actual_resolution def test_function_validate_layer_instance(): - is_valid, except_str = validate_layer_instance(Albedo()) - assert is_valid is False is_valid, except_str = validate_layer_instance('t') assert is_valid is False - is_valid, except_str = validate_layer_instance('Layer()') + is_valid, except_str = validate_layer_instance(EsaWorldCoverClass.BUILT_UP) assert is_valid is False - is_valid, except_str = validate_layer_instance('OpenStreetMap()') + is_valid, except_str = validate_layer_instance(OpenStreetMap()) assert is_valid is False - is_valid, except_str = validate_layer_instance('Albedo(spatial_resolution = 2)') + is_valid, except_str = validate_layer_instance(Albedo(spatial_resolution = 2)) assert is_valid is False -def validate_layer_instance(obj_string): +def validate_layer_instance(obj): is_valid = True except_str = None - obj_eval = None - - if not type(obj_string) == str: - is_valid = False - except_str = "Specified object '%s' must be specified as a string." % obj_string - return is_valid, except_str - try: - obj_eval = eval(obj_string) - except: + if not obj.__class__.__bases__[0] == Layer: is_valid = False - except_str = "Specified object '%s' is not a class instance." % obj_string - return is_valid, except_str - - if not type(obj_eval).__bases__[0] == Layer: - is_valid = False - except_str = "Specified object '%s' is not a valid Layer class instance." % obj_string - elif not hasattr(obj_eval, 'spatial_resolution'): - is_valid = False - except_str = "Specified class '%s' does not have the spatial_resolution property." % obj_string - elif not obj_string.find('spatial_resolution') == -1: - is_valid = False - except_str = "Do not specify spatial_resolution property value in object '%s'." % obj_string - elif obj_eval.spatial_resolution is None: - is_valid = False - except_str = "Class signature cannot specify None for default value for class." + except_str = "Specified object '%s' is not a valid Layer class instance." % obj + else: + cls = get_class_from_instance(obj) + cls_name = type(cls).__name__ + if not hasattr(obj, 'spatial_resolution'): + is_valid = False + except_str = "Class '%s' does not have spatial_resolution property." % cls_name + elif not obj.spatial_resolution == cls.spatial_resolution: + is_valid = False + except_str = "Do not specify spatial_resolution property value for class '%s'." % cls_name + elif cls.spatial_resolution is None: + is_valid = False + except_str = "Signature of class %s must specify a non-null default value for spatial_resolution. Please correct." % cls_name return is_valid, except_str -def evaluate_layer(layer, expected_resolution): - data = layer.get_data(BBOX) - actual_estimated_resolution = get_spatial_resolution_estimate(data) - assert expected_resolution == actual_estimated_resolution +def get_class_from_instance(obj): + cls = obj.__class__() + return cls -def get_spatial_resolution_estimate(data): +def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) + y_min = data.coords['y'].values.min() + y_max = data.coords['y'].values.max() + + crs_string = data.rio.crs.data['init'] + crs = CRS.from_string(crs_string) + crs_unit = crs.axis_info[0].unit_name - y_min = data['y'].values.min() - y_max = data['y'].values.max() - - crs = CRS.from_string(data.crs) - crs_units = crs.axis_info[0].unit_name - if crs_units == 'metre': - y_diff = y_max - y_min - elif crs_units == 'foot': - feet_to_meter = 0.3048 - y_diff = (y_max - y_min) * feet_to_meter - elif crs_units == 'degree': - lat1 = y_min - lat2 = y_max - lon1 = data['x'].values.min() - lon2 = lon1 - y_diff = get_distance_between_geocoordinates(lat1, lon1, lat2, lon2) + if crs_unit == 'metre': + diff_distance = y_max - y_min else: - raise Exception('Unhandled projection units: %s' % crs_units) + raise Exception('Unhandled projection units: %s for projection: %s' % (crs_unit, crs_string)) - ry = round(y_diff / y_cells) + estimated_actual_resolution = round(diff_distance / y_cells) - return ry + return estimated_actual_resolution diff --git a/tests/tools/spatial_tools.py b/tests/tools/spatial_tools.py deleted file mode 100644 index 73ab093..0000000 --- a/tests/tools/spatial_tools.py +++ /dev/null @@ -1,16 +0,0 @@ -import math - -EARTH_RADIUS = 6378.137 # Radius of earth in KM - -def get_distance_between_geocoordinates(lat1, lon1, lat2, lon2): - dLat = lat2 * math.pi / 180 - lat1 * math.pi / 180 - dLon = lon2 * math.pi / 180 - lon1 * math.pi / 180 - a = math.sin(dLat/2) * math.sin(dLat/2) + \ - math.cos(lat1 * math.pi / 180) * \ - math.cos(lat2 * math.pi / 180) * \ - math.sin(dLon/2) * math.sin(dLon/2) - c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) - d = EARTH_RADIUS * c - return d * 1000 # meters - - From 8e88425d5d988d49aeba37f8d7071649239f7942 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Wed, 4 Sep 2024 09:52:49 -0700 Subject: [PATCH 28/32] Cleaned up and expanded resolution testing --- .../layers/high_land_surface_temperature.py | 5 +- city_metrix/layers/layer.py | 4 +- .../conftest.py | 12 ++- .../test_write_layers_to_qgis_files.py | 80 +++++++++++------- tests/test_layer_parameters.py | 84 ++++++++++++------- tests/tools/general_tools.py | 11 ++- 6 files changed, 127 insertions(+), 69 deletions(-) diff --git a/city_metrix/layers/high_land_surface_temperature.py b/city_metrix/layers/high_land_surface_temperature.py index d3943e9..440eb02 100644 --- a/city_metrix/layers/high_land_surface_temperature.py +++ b/city_metrix/layers/high_land_surface_temperature.py @@ -10,17 +10,18 @@ class HighLandSurfaceTemperature(Layer): THRESHOLD_ADD = 3 - def __init__(self, start_date="2013-01-01", end_date="2023-01-01", **kwargs): + def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resolution=30, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date + self.spatial_resolution = spatial_resolution def get_data(self, bbox): hottest_date = self.get_hottest_date(bbox) start_date = (hottest_date - datetime.timedelta(days=45)).strftime("%Y-%m-%d") end_date = (hottest_date + datetime.timedelta(days=45)).strftime("%Y-%m-%d") - lst = LandSurfaceTemperature(start_date, end_date).get_data(bbox) + lst = LandSurfaceTemperature(start_date, end_date, self.spatial_resolution).get_data(bbox) lst_mean = lst.mean(dim=['x', 'y']) high_lst = lst.where(lst >= (lst_mean + self.THRESHOLD_ADD)) return high_lst diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 86590c8..c86e9fe 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -274,10 +274,12 @@ def get_image_collection( :param name: optional name to print while reporting progress :return: """ + if scale is None: + raise Exception("Spatial_resolution cannot be None.") crs = get_utm_zone_epsg(bbox) - # See link regarding bug in crs specification shttps://github.com/google/Xee/issues/118 + # See link regarding bug in crs specification https://github.com/google/Xee/issues/118 ds = xr.open_dataset( image_collection, engine='ee', diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py index f8727f5..7ca4581 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py +++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py @@ -12,13 +12,17 @@ # Values should normally be set to False in order to avoid unnecessary execution. RUN_DUMPS = False -# Specify None to write to a temporary default folder otherwise specify a valid custom target path. -CUSTOM_DUMP_DIRECTORY = None +# Multiplier applied to the default spatial_resolution of the layer +# Use value of 1 for default resolution. +SPATIAL_RESOLUTION_MULTIPLIER = 1 # Both the tests and QGIS file are implemented for the same bounding box in Brazil. COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 +# Specify None to write to a temporary default folder otherwise specify a valid custom target path. +CUSTOM_DUMP_DIRECTORY = None + def pytest_configure(config): qgis_project_file = 'layers_for_br_lauro_de_freitas.qgz' @@ -36,6 +40,10 @@ def pytest_configure(config): def target_folder(): return get_target_folder_path() +@pytest.fixture +def target_spatial_resolution_multiplier(): + return SPATIAL_RESOLUTION_MULTIPLIER + @pytest.fixture def bbox_info(): bbox = namedtuple('bbox', ['bounds', 'country']) diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py index 5e0efb9..5e89328 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py +++ b/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py @@ -24,122 +24,138 @@ WorldPop, Layer ) from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated +from ...tools.general_tools import get_class_default_spatial_resolution + @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_albedo(target_folder, bbox_info): +def test_write_albedo(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'albedo.tif') - Albedo().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(Albedo()) + Albedo(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_alos_dsm(target_folder, bbox_info): +def test_write_alos_dsm(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'alos_dsm.tif') - AlosDSM().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(AlosDSM()) + AlosDSM(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_average_net_building_height(target_folder, bbox_info): +def test_write_average_net_building_height(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'average_net_building_height.tif') - AverageNetBuildingHeight().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(AverageNetBuildingHeight()) + AverageNetBuildingHeight(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_esa_world_cover(target_folder, bbox_info): +def test_write_esa_world_cover(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'esa_world_cover.tif') - EsaWorldCover().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(EsaWorldCover()) + EsaWorldCover(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_high_land_surface_temperature(target_folder, bbox_info): +def test_write_high_land_surface_temperature(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'high_land_surface_temperature.tif') - HighLandSurfaceTemperature().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(HighLandSurfaceTemperature()) + HighLandSurfaceTemperature(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_land_surface_temperature(target_folder, bbox_info): +def test_write_land_surface_temperature(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'land_surface_temperature.tif') - LandSurfaceTemperature().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(LandSurfaceTemperature()) + LandSurfaceTemperature(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) # TODO Class is no longer used, but may be useful later # @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -# def test_write_landsat_collection_2(target_folder, bbox_info): +# def test_write_landsat_collection_2(target_folder, bbox_info, target_spatial_resolution_multiplier): # file_path = prep_output_path(target_folder, 'landsat_collection2.tif') # bands = ['green'] # LandsatCollection2(bands).write(bbox_info.bounds, file_path, tile_degrees=None) # assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_nasa_dem(target_folder, bbox_info): +def test_write_nasa_dem(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'nasa_dem.tif') - NasaDEM().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(NasaDEM()) + NasaDEM(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_natural_areas(target_folder, bbox_info): +def test_write_natural_areas(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'natural_areas.tif') - NaturalAreas().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(NaturalAreas()) + NaturalAreas(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_ndvi_sentinel2_gee(target_folder, bbox_info): +def test_write_ndvi_sentinel2_gee(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'ndvi_sentinel2_gee.tif') - NdviSentinel2(year=2023).write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(NdviSentinel2()) + NdviSentinel2(year=2023, spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_openbuildings(target_folder, bbox_info): +def test_write_openbuildings(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'open_buildings.tif') OpenBuildings(bbox_info.country).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) # TODO Class write is not functional. Is class still needed or have we switched to overture? # @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -# def test_write_open_street_map(target_folder, bbox_info): +# def test_write_open_street_map(target_folder, bbox_info, target_spatial_resolution_multiplier): # file_path = prep_output_path(target_folder, 'open_street_map.tif') # OpenStreetMap().write(bbox_info.bounds, file_path, tile_degrees=None) # assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_overture_buildings(target_folder, bbox_info): +def test_write_overture_buildings(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'overture_buildings.tif') OvertureBuildings().write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) # TODO Class is no longer used, but may be useful later # @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -# def test_write_sentinel_2_level2(target_folder, bbox_info): +# def test_write_sentinel_2_level2(target_folder, bbox_info, target_spatial_resolution_multiplier): # file_path = prep_output_path(target_folder, 'sentinel_2_level2.tif') # sentinel_2_bands = ["green"] # Sentinel2Level2(sentinel_2_bands).write(bbox_info.bounds, file_path, tile_degrees=None) # assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_smart_surface_lulc(target_folder, bbox_info): +def test_write_smart_surface_lulc(target_folder, bbox_info, target_spatial_resolution_multiplier): + # Note: spatial_resolution not implemented for this raster class file_path = prep_output_path(target_folder, 'smart_surface_lulc.tif') SmartSurfaceLULC().write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_tree_canopy_height(target_folder, bbox_info): +def test_write_tree_canopy_height(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'tree_canopy_height.tif') - TreeCanopyHeight().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(TreeCanopyHeight()) + TreeCanopyHeight(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_tree_cover(target_folder, bbox_info): +def test_write_tree_cover(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'tree_cover.tif') - TreeCover().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(TreeCover()) + TreeCover(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_urban_land_use(target_folder, bbox_info): +def test_write_urban_land_use(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'urban_land_use.tif') - UrbanLandUse().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(UrbanLandUse()) + UrbanLandUse(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_world_pop(target_folder, bbox_info): +def test_write_world_pop(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'world_pop.tif') - WorldPop().write(bbox_info.bounds, file_path, tile_degrees=None) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(WorldPop()) + WorldPop(spatial_resolution=target_resolution).write(bbox_info.bounds, file_path, tile_degrees=None) assert verify_file_is_populated(file_path) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 2238bad..171b338 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -14,9 +14,10 @@ TreeCanopyHeight, TreeCover, UrbanLandUse, - WorldPop, OpenStreetMap + WorldPop, OpenStreetMap, HighLandSurfaceTemperature ) from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 +from tests.tools.general_tools import get_class_from_instance, get_class_default_spatial_resolution """ Evaluation of spatial_resolution property @@ -33,86 +34,90 @@ def test_albedo_spatial_resolution(): class_instance = Albedo() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_alos_dsm_spatial_resolution(): class_instance = AlosDSM() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_average_net_building_height_spatial_resolution(): class_instance = AverageNetBuildingHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_built_up_height_spatial_resolution(): class_instance = BuiltUpHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_esa_world_cover_spatial_resolution(): class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) + assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + +def test_high_land_surface_temperature_spatial_resolution(): + class_instance = HighLandSurfaceTemperature() + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_land_surface_temperature_spatial_resolution(): class_instance = LandSurfaceTemperature() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_nasa_dem_spatial_resolution(): class_instance = NasaDEM() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_natural_areas_spatial_resolution(): class_instance = NaturalAreas() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_ndvi_sentinel2_spatial_resolution(): class_instance = NdviSentinel2(year=2023) - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_tree_canopy_height_spatial_resolution(): class_instance = TreeCanopyHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_tree_cover_spatial_resolution(): class_instance = TreeCover() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_urban_land_use_spatial_resolution(): class_instance = UrbanLandUse() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution def test_world_pop_spatial_resolution(): class_instance = WorldPop() - doubled_default_resolution, actual_estimated_resolution = evaluate_resolution__property(class_instance) + doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution -def evaluate_resolution__property(obj): - is_valid, except_str = validate_layer_instance(obj) - if is_valid is False: - raise Exception(except_str) - - # Double the default scale for testing - cls = get_class_from_instance(obj) - doubled_default_resolution = 2 * cls.spatial_resolution - obj.spatial_resolution=doubled_default_resolution +def test_halved_spatial_resolution(): + class_instance = Albedo() + halved_default_resolution = get_class_default_spatial_resolution(class_instance) / 2 + class_instance.spatial_resolution=halved_default_resolution - data = obj.get_data(BBOX) + expected_resolution = halved_default_resolution + estimated_actual_resolution = get_modified_resolution_data(class_instance) - expected_resolution = doubled_default_resolution - estimated_actual_resolution = estimate_spatial_resolution(data) + assert expected_resolution == estimated_actual_resolution - return expected_resolution, estimated_actual_resolution +def test_null_spatial_resolution(): + class_instance = Albedo() + class_instance.spatial_resolution=None + with pytest.raises(Exception) as e_info: + get_modified_resolution_data(class_instance) def test_function_validate_layer_instance(): is_valid, except_str = validate_layer_instance('t') @@ -124,6 +129,27 @@ def test_function_validate_layer_instance(): is_valid, except_str = validate_layer_instance(Albedo(spatial_resolution = 2)) assert is_valid is False +def evaluate_doubled_resolution_property(obj): + is_valid, except_str = validate_layer_instance(obj) + if is_valid is False: + raise Exception(except_str) + + # Double the default scale for testing + doubled_default_resolution = 2 * get_class_default_spatial_resolution(obj) + obj.spatial_resolution=doubled_default_resolution + + expected_resolution = doubled_default_resolution + estimated_actual_resolution = get_modified_resolution_data(obj) + + return expected_resolution, estimated_actual_resolution + +def get_modified_resolution_data(obj): + data = obj.get_data(BBOX) + estimated_actual_resolution = estimate_spatial_resolution(data) + + return estimated_actual_resolution + + def validate_layer_instance(obj): is_valid = True except_str = None @@ -146,10 +172,6 @@ def validate_layer_instance(obj): return is_valid, except_str -def get_class_from_instance(obj): - cls = obj.__class__() - return cls - def estimate_spatial_resolution(data): y_cells = float(data['y'].size - 1) y_min = data.coords['y'].values.min() diff --git a/tests/tools/general_tools.py b/tests/tools/general_tools.py index 99bd894..d3e56b8 100644 --- a/tests/tools/general_tools.py +++ b/tests/tools/general_tools.py @@ -54,4 +54,13 @@ def convert_ratio_to_percentage(data): source_crs = data.rio.crs values_as_percent.rio.write_crs(source_crs, inplace=True) - return values_as_percent \ No newline at end of file + return values_as_percent + +def get_class_default_spatial_resolution(obj): + cls = get_class_from_instance(obj) + default_spatial_resolution = cls.spatial_resolution + return default_spatial_resolution + +def get_class_from_instance(obj): + cls = obj.__class__() + return cls \ No newline at end of file From 8ff3340bab706fffd93a25812c3e1c20e225350a Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Wed, 4 Sep 2024 16:28:51 -0700 Subject: [PATCH 29/32] Added docstring --- city_metrix/layers/high_land_surface_temperature.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/city_metrix/layers/high_land_surface_temperature.py b/city_metrix/layers/high_land_surface_temperature.py index 440eb02..5faa2ae 100644 --- a/city_metrix/layers/high_land_surface_temperature.py +++ b/city_metrix/layers/high_land_surface_temperature.py @@ -8,6 +8,12 @@ class HighLandSurfaceTemperature(Layer): + """ + Attributes: + start_date: starting date for data retrieval + end_date: ending date for data retrieval + spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + """ THRESHOLD_ADD = 3 def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resolution=30, **kwargs): From efeed2e99237c899e03e345294dbd90bf91a818a Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 5 Sep 2024 14:18:26 -0700 Subject: [PATCH 30/32] Added test of downsized image values --- city_metrix/layers/ndvi_sentinel2_gee.py | 7 +- tests/test_layer_parameters.py | 220 ++++++++++++++++------- 2 files changed, 160 insertions(+), 67 deletions(-) diff --git a/city_metrix/layers/ndvi_sentinel2_gee.py b/city_metrix/layers/ndvi_sentinel2_gee.py index 89b4c09..e49b3d3 100644 --- a/city_metrix/layers/ndvi_sentinel2_gee.py +++ b/city_metrix/layers/ndvi_sentinel2_gee.py @@ -42,8 +42,7 @@ def calculate_ndvi(image): ndvi_mosaic = ndvi.qualityMosaic('NDVI') ic = ee.ImageCollection(ndvi_mosaic) - ndvi_data = get_image_collection(ic, bbox, self.spatial_resolution, "NDVI") + ndvi_data = (get_image_collection(ic, bbox, self.spatial_resolution, "NDVI") + .NDVI) - xdata = ndvi_data.to_dataarray() - - return xdata + return ndvi_data diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 171b338..e085c73 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,4 +1,6 @@ import pytest +import numpy as np +from skimage.metrics import structural_similarity as ssim from pyproj import CRS from city_metrix.layers import ( Layer, @@ -30,94 +32,137 @@ COUNTRY_CODE_FOR_BBOX = 'BRA' BBOX = BBOX_BRA_LAURO_DE_FREITAS_1 -RESOLUTION_TOLERANCE = 1 +RESOLUTION_COMPARISON_TOLERANCE = 1 +DOWNSIZE_FACTOR = 2 -def test_albedo_spatial_resolution(): +def test_albedo_downsampled_spatial_resolution(): class_instance = Albedo() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) -def test_alos_dsm_spatial_resolution(): + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances + +def test_alos_dsm_downsampled_spatial_resolution(): class_instance = AlosDSM() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) + + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances -def test_average_net_building_height_spatial_resolution(): +def test_average_net_building_height_downsampled_spatial_resolution(): class_instance = AverageNetBuildingHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) -def test_built_up_height_spatial_resolution(): + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances + +def test_built_up_height_downsampled_spatial_resolution(): class_instance = BuiltUpHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) + + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances -def test_esa_world_cover_spatial_resolution(): +def test_esa_world_cover_downsampled_spatial_resolution(): class_instance = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP) - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) -def test_high_land_surface_temperature_spatial_resolution(): + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances + +def test_high_land_surface_temperature_downsampled_spatial_resolution(): class_instance = HighLandSurfaceTemperature() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) + + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances -def test_land_surface_temperature_spatial_resolution(): +def test_land_surface_temperature_downsampled_spatial_resolution(): class_instance = LandSurfaceTemperature() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) -def test_nasa_dem_spatial_resolution(): + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances + +def test_nasa_dem_downsampled_spatial_resolution(): class_instance = NasaDEM() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) + + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances -def test_natural_areas_spatial_resolution(): +def test_natural_areas_downsampled_spatial_resolution(): class_instance = NaturalAreas() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) -def test_ndvi_sentinel2_spatial_resolution(): + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances + +def test_ndvi_sentinel2_downsampled_spatial_resolution(): class_instance = NdviSentinel2(year=2023) - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) + + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances -def test_tree_canopy_height_spatial_resolution(): +def test_tree_canopy_height_downsampled_spatial_resolution(): class_instance = TreeCanopyHeight() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) -def test_tree_cover_spatial_resolution(): + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances + +def test_tree_cover_downsampled_spatial_resolution(): class_instance = TreeCover() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) + + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances -def test_urban_land_use_spatial_resolution(): +def test_urban_land_use_downsampled_spatial_resolution(): class_instance = UrbanLandUse() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) -def test_world_pop_spatial_resolution(): + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances + +def test_world_pop_downsampled_spatial_resolution(): class_instance = WorldPop() - doubled_default_resolution, actual_estimated_resolution = evaluate_doubled_resolution_property(class_instance) - assert pytest.approx(doubled_default_resolution, rel=RESOLUTION_TOLERANCE) == actual_estimated_resolution + default_res_data, downsized_res_data, target_downsized_res, estimated_actual_res = get_sample_data(class_instance) + downsizing_is_within_tolerances = evaluate_raster_value(default_res_data, downsized_res_data) + + assert pytest.approx(target_downsized_res, rel=RESOLUTION_COMPARISON_TOLERANCE) == estimated_actual_res + assert downsizing_is_within_tolerances -def test_halved_spatial_resolution(): +def test_halved_up_sampled_spatial_resolution(): class_instance = Albedo() halved_default_resolution = get_class_default_spatial_resolution(class_instance) / 2 - class_instance.spatial_resolution=halved_default_resolution expected_resolution = halved_default_resolution - estimated_actual_resolution = get_modified_resolution_data(class_instance) + modified_data = get_modified_resolution_data(class_instance, halved_default_resolution) + estimated_actual_resolution = estimate_spatial_resolution(modified_data) assert expected_resolution == estimated_actual_resolution def test_null_spatial_resolution(): class_instance = Albedo() - class_instance.spatial_resolution=None + spatial_resolution=None with pytest.raises(Exception) as e_info: - get_modified_resolution_data(class_instance) + get_modified_resolution_data(class_instance, spatial_resolution) def test_function_validate_layer_instance(): is_valid, except_str = validate_layer_instance('t') @@ -129,26 +174,31 @@ def test_function_validate_layer_instance(): is_valid, except_str = validate_layer_instance(Albedo(spatial_resolution = 2)) assert is_valid is False -def evaluate_doubled_resolution_property(obj): - is_valid, except_str = validate_layer_instance(obj) +def get_sample_data(class_instance): + is_valid, except_str = validate_layer_instance(class_instance) if is_valid is False: raise Exception(except_str) - # Double the default scale for testing - doubled_default_resolution = 2 * get_class_default_spatial_resolution(obj) - obj.spatial_resolution=doubled_default_resolution + default_res = get_class_default_spatial_resolution(class_instance) + downsized_resolution = DOWNSIZE_FACTOR * default_res - expected_resolution = doubled_default_resolution - estimated_actual_resolution = get_modified_resolution_data(obj) + downsized_res_data = get_modified_resolution_data(class_instance, downsized_resolution) + default_res_data = get_modified_resolution_data(class_instance, default_res) - return expected_resolution, estimated_actual_resolution + estimated_actual_resolution = estimate_spatial_resolution(downsized_res_data) -def get_modified_resolution_data(obj): - data = obj.get_data(BBOX) - estimated_actual_resolution = estimate_spatial_resolution(data) + return default_res_data, downsized_res_data, downsized_resolution, estimated_actual_resolution + +def get_crs_from_image_data(image_data): + crs_string = image_data.rio.crs.data['init'] + crs = CRS.from_string(crs_string) + return crs - return estimated_actual_resolution +def get_modified_resolution_data(class_instance, spatial_resolution): + class_instance.spatial_resolution = spatial_resolution + data = class_instance.get_data(BBOX) + return data def validate_layer_instance(obj): is_valid = True @@ -177,15 +227,59 @@ def estimate_spatial_resolution(data): y_min = data.coords['y'].values.min() y_max = data.coords['y'].values.max() - crs_string = data.rio.crs.data['init'] - crs = CRS.from_string(crs_string) + crs = get_crs_from_image_data(data) crs_unit = crs.axis_info[0].unit_name if crs_unit == 'metre': diff_distance = y_max - y_min else: - raise Exception('Unhandled projection units: %s for projection: %s' % (crs_unit, crs_string)) + raise Exception('Unhandled projection units: %s for projection: %s' % (crs_unit, crs.srs)) estimated_actual_resolution = round(diff_distance / y_cells) return estimated_actual_resolution + +def get_populate_ratio(dataset): + raw_data_size = dataset.values.size + populated_raw_data = dataset.values[(~np.isnan(dataset.values)) & (dataset.values > 0)] + populated_data_raw_size = populated_raw_data.size + populated_raw_data_ratio = populated_data_raw_size/raw_data_size + return populated_raw_data_ratio + +def evaluate_raster_value(raw_data, downsized_data): + # Below values where determined through trial and error evaluation of results in QGIS + ratio_tolerance = 0.2 + normalized_rmse_tolerance = 0.3 + ssim_index_tolerance = 0.6 + + populated_raw_data_ratio = get_populate_ratio(raw_data) + populated_downsized_data_ratio = get_populate_ratio(raw_data) + ratio_eval = are_numbers_within_tolerance(populated_raw_data_ratio, populated_downsized_data_ratio, ratio_tolerance) + + filled_raw_data = raw_data.fillna(0) + filled_downsized_data = downsized_data.fillna(0) + + # Resample raw_data to match the resolution of downsized_data + resampled_raw_data = filled_raw_data.interp_like(filled_downsized_data).fillna(0) + + # Convert xarray DataArrays to numpy arrays + processed_downsized_data_np = filled_downsized_data.values + processed_raw_data_np = resampled_raw_data.values + + # Calculate and evaluate normalized Mean Squared Error (MSE) + max_val = processed_downsized_data_np.max() \ + if processed_downsized_data_np.max() > processed_raw_data_np.max() else processed_raw_data_np.max() + normalized_rmse = np.sqrt(np.mean(np.square(processed_downsized_data_np - processed_raw_data_np))) / max_val + matching_rmse = True if normalized_rmse < normalized_rmse_tolerance else False + + # Calculate and evaluate Structural Similarity Index (SSIM) + ssim_index, _ = ssim(processed_downsized_data_np, processed_raw_data_np, full=True, data_range=max_val) + matching_ssim = True if ssim_index > ssim_index_tolerance else False + + results_match = True if (ratio_eval & matching_rmse & matching_ssim) else False + + return results_match + +def are_numbers_within_tolerance(num1, num2, tolerance): + diff = abs(num1 - num2) + return True if diff <= tolerance else False From 809c39d6f308d38cd7deb1895dede341c1015fb2 Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 5 Sep 2024 14:26:51 -0700 Subject: [PATCH 31/32] Added missing packages --- .github/requirements.txt | 1 + environment.yml | 1 + setup.py | 1 + 3 files changed, 3 insertions(+) diff --git a/.github/requirements.txt b/.github/requirements.txt index 15d8543..90f10fe 100644 --- a/.github/requirements.txt +++ b/.github/requirements.txt @@ -16,5 +16,6 @@ geemap==0.32.0 pip==23.3.1 boto3==1.34.124 scikit-learn==1.5.0 +scikit-image==0.24.0 overturemaps==0.6.0 exactextract==0.2.0.dev252 diff --git a/environment.yml b/environment.yml index 2434184..18ac340 100644 --- a/environment.yml +++ b/environment.yml @@ -22,6 +22,7 @@ dependencies: - pip=23.3.1 - boto3=1.34.124 - scikit-learn=1.5.0 + - scikit-image==0.24.0 - exactextract=0.2.0.dev252 - pip: - overturemaps==0.6.0 diff --git a/setup.py b/setup.py index 179b595..ff88cbd 100644 --- a/setup.py +++ b/setup.py @@ -31,5 +31,6 @@ "exactextract<=0.2.0.dev252", "overturemaps", "scikit-learn>=1.5.0", + "scikit-image>=0.24.0" ], ) From ee0d8d171804bcb52dfc72c0cf34952bc381b6eb Mon Sep 17 00:00:00 2001 From: Kenn Cartier Date: Thu, 5 Sep 2024 17:41:58 -0700 Subject: [PATCH 32/32] Added test documentation and switch to quadratic interpolation for image resampling --- tests/test_layer_parameters.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index e085c73..90942f7 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -2,6 +2,8 @@ import numpy as np from skimage.metrics import structural_similarity as ssim from pyproj import CRS +from xrspatial import quantile + from city_metrix.layers import ( Layer, Albedo, @@ -254,19 +256,28 @@ def evaluate_raster_value(raw_data, downsized_data): populated_raw_data_ratio = get_populate_ratio(raw_data) populated_downsized_data_ratio = get_populate_ratio(raw_data) - ratio_eval = are_numbers_within_tolerance(populated_raw_data_ratio, populated_downsized_data_ratio, ratio_tolerance) + diff = abs(populated_raw_data_ratio - populated_downsized_data_ratio) + ratio_eval = True if diff <= ratio_tolerance else False filled_raw_data = raw_data.fillna(0) filled_downsized_data = downsized_data.fillna(0) - # Resample raw_data to match the resolution of downsized_data - resampled_raw_data = filled_raw_data.interp_like(filled_downsized_data).fillna(0) + # Resample raw_data to match the resolution of the downsized_data. + # This operation is necessary in order to use RMSE and SSIM since they require matching array dimensions. + # Interpolation is set to quadratic method to intentionally not match method used by xee. (TODO Find documentation) + # For future releases of this method, we may need to control the interpolation to match what was used + # for a specific layer. + # Note: Initial investigations using the unresampled-raw and downsized data with evaluation by + # mean value, quantiles,and standard deviation were unsuccessful due to false failures on valid downsized images. + resampled_filled_raw_data = (filled_raw_data + .interp_like(filled_downsized_data, method='quadratic') + .fillna(0)) # Convert xarray DataArrays to numpy arrays + processed_raw_data_np = resampled_filled_raw_data.values processed_downsized_data_np = filled_downsized_data.values - processed_raw_data_np = resampled_raw_data.values - # Calculate and evaluate normalized Mean Squared Error (MSE) + # Calculate and evaluate normalized Root Mean Squared Error (RMSE) max_val = processed_downsized_data_np.max() \ if processed_downsized_data_np.max() > processed_raw_data_np.max() else processed_raw_data_np.max() normalized_rmse = np.sqrt(np.mean(np.square(processed_downsized_data_np - processed_raw_data_np))) / max_val @@ -279,7 +290,3 @@ def evaluate_raster_value(raw_data, downsized_data): results_match = True if (ratio_eval & matching_rmse & matching_ssim) else False return results_match - -def are_numbers_within_tolerance(num1, num2, tolerance): - diff = abs(num1 - num2) - return True if diff <= tolerance else False