diff --git a/README.md b/README.md index 56ae15e6..4af03c1a 100644 --- a/README.md +++ b/README.md @@ -43,7 +43,7 @@ To run the module, 1. You need access to Google Earth Engine 2. Install - 3. If you want to use the ERA5 layer, you need to install the [Climate Data Store (CDS) Application Program Interface (API)](https://cds.climate.copernicus.eu/api-how-to) + 3. If you want to use the ERA5 layer, you need to install the [Climate Data Store (CDS) Application Program Interface (API)](https://cds.climate.copernicus.eu/how-to-api) ### Interactive development diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py index 341786b2..5cb7aaaf 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -2,7 +2,8 @@ import xarray from dask.diagnostics import ProgressBar -from .layer import Layer, get_utm_zone_epsg, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection, set_resampling_for_continuous_raster + class Albedo(Layer): """ @@ -10,17 +11,20 @@ class Albedo(Layer): 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) + resampling_method: interpolation method used by Google Earth Engine. Default is 'bilinear'. All options are: ('bilinear', 'bicubic', None). 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): + def __init__(self, start_date="2021-01-01", end_date="2022-01-01", spatial_resolution:int=10, + resampling_method:str='bilinear', threshold=None, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date self.spatial_resolution = spatial_resolution + self.resampling_method = resampling_method self.threshold = threshold - def get_data(self, bbox): + def get_data(self, bbox: tuple[float, float, float, float]): S2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") S2C = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY") @@ -34,19 +38,33 @@ def get_data(self, bbox): def mask_and_count_clouds(s2wc, geom): s2wc = ee.Image(s2wc) geom = ee.Geometry(geom.geometry()) - is_cloud = ee.Image(s2wc.get('cloud_mask')).gt(MAX_CLOUD_PROB).rename('is_cloud') + is_cloud = (ee.Image(s2wc.get('cloud_mask')) + .gt(MAX_CLOUD_PROB) + .rename('is_cloud') + ) + nb_cloudy_pixels = is_cloud.reduceRegion( reducer=ee.Reducer.sum().unweighted(), geometry=geom, scale=self.spatial_resolution, maxPixels=1e9 ) - return s2wc.updateMask(is_cloud.eq(0)).set('nb_cloudy_pixels', - nb_cloudy_pixels.getNumber('is_cloud')).divide(10000) + mask = (s2wc + .updateMask(is_cloud.eq(0)) + .set('nb_cloudy_pixels',nb_cloudy_pixels.getNumber('is_cloud')) + .divide(10000) + ) + + return mask def mask_clouds_and_rescale(im): - clouds = ee.Image(im.get('cloud_mask')).select('probability') - return im.updateMask(clouds.lt(MAX_CLOUD_PROB)).divide(10000) + clouds = ee.Image(im.get('cloud_mask') + ).select('probability') + mask = im.updateMask(clouds + .lt(MAX_CLOUD_PROB) + ).divide(10000) + + return mask def get_masked_s2_collection(roi, start, end): criteria = (ee.Filter.And( @@ -55,23 +73,29 @@ def get_masked_s2_collection(roi, start, end): )) s2 = S2.filter(criteria) # .select('B2','B3','B4','B8','B11','B12') s2c = S2C.filter(criteria) - s2_with_clouds = (ee.Join.saveFirst('cloud_mask').apply(**{ - 'primary': ee.ImageCollection(s2), - 'secondary': ee.ImageCollection(s2c), - 'condition': ee.Filter.equals(**{'leftField': 'system:index', 'rightField': 'system:index'}) - })) + s2_with_clouds = ( + ee.Join.saveFirst('cloud_mask') + .apply(**{ + 'primary': ee.ImageCollection(s2), + 'secondary': ee.ImageCollection(s2c), + 'condition': ee.Filter.equals(**{'leftField': 'system:index', 'rightField': 'system:index'}) + }) + ) def _mcc(im): return mask_and_count_clouds(im, roi) # s2_with_clouds=ee.ImageCollection(s2_with_clouds).map(_mcc) # s2_with_clouds=s2_with_clouds.limit(image_limit,'nb_cloudy_pixels') - s2_with_clouds = ee.ImageCollection(s2_with_clouds).map( - mask_clouds_and_rescale) # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE') - return ee.ImageCollection(s2_with_clouds) + s2_with_clouds = (ee.ImageCollection(s2_with_clouds) + .map(mask_clouds_and_rescale) + ) # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE') - # calculate albedo for images + s2_with_clouds_ic = ee.ImageCollection(s2_with_clouds) + return s2_with_clouds_ic + + # calculate albedo for images # weights derived from # S. Bonafoni and A. Sekertekin, "Albedo Retrieval From Sentinel-2 by New Narrow-to-Broadband Conversion Coefficients," in IEEE Geoscience and Remote Sensing Letters, vol. 17, no. 9, pp. 1618-1622, Sept. 2020, doi: 10.1109/LGRS.2020.2967085. def calc_s2_albedo(image): @@ -89,20 +113,35 @@ def calc_s2_albedo(image): 'SWIR1': image.select('B11'), 'SWIR2': image.select('B12') } - return image.expression(S2_ALBEDO_EQN, config).double().rename('albedo') + + albedo = image.expression(S2_ALBEDO_EQN, config).double().rename('albedo') + + return albedo ## S2 MOSAIC AND ALBEDO dataset = get_masked_s2_collection(ee.Geometry.BBox(*bbox), self.start_date, self.end_date) s2_albedo = dataset.map(calc_s2_albedo) - albedo_mean = s2_albedo.reduce(ee.Reducer.mean()) - data = (get_image_collection( - ee.ImageCollection(albedo_mean), bbox, self.spatial_resolution, "albedo") - .albedo_mean) + albedo_mean = (s2_albedo + .map(lambda x: + set_resampling_for_continuous_raster(x, + self.resampling_method, + self.spatial_resolution, + bbox + ) + ) + .reduce(ee.Reducer.mean()) + ) + + albedo_mean_ic = ee.ImageCollection(albedo_mean) + data = get_image_collection( + albedo_mean_ic, + bbox, + self.spatial_resolution, + "albedo" + ).albedo_mean if self.threshold is not None: return data.where(data < self.threshold) return data - - diff --git a/city_metrix/layers/alos_dsm.py b/city_metrix/layers/alos_dsm.py index 70000eb5..f9601eb8 100644 --- a/city_metrix/layers/alos_dsm.py +++ b/city_metrix/layers/alos_dsm.py @@ -2,26 +2,45 @@ import xee import xarray as xr -from .layer import Layer, get_image_collection + +from .layer import Layer, get_image_collection, set_resampling_for_continuous_raster class AlosDSM(Layer): """ Attributes: spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + resampling_method: interpolation method used by Google Earth Engine. Default is 'bilinear'. All options are: ('bilinear', 'bicubic', None). """ - def __init__(self, spatial_resolution=30, **kwargs): + def __init__(self, spatial_resolution:int=30, resampling_method:str='bilinear', **kwargs): super().__init__(**kwargs) self.spatial_resolution = spatial_resolution + self.resampling_method = resampling_method + + def get_data(self, bbox: tuple[float, float, float, float]): + alos_dsm = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2") + + alos_dsm_ic = ee.ImageCollection( + alos_dsm + .filterBounds(ee.Geometry.BBox(*bbox)) + .select('DSM') + .map(lambda x: + set_resampling_for_continuous_raster(x, + self.resampling_method, + self.spatial_resolution, + bbox + ) + ) + .mean() + ) + - def get_data(self, bbox): - dataset = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2") - alos_dsm = ee.ImageCollection(dataset - .filterBounds(ee.Geometry.BBox(*bbox)) - .select('DSM') - .mean() - ) - data = get_image_collection(alos_dsm, bbox, self.spatial_resolution, "ALOS DSM").DSM + data = get_image_collection( + alos_dsm_ic, + 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 d0f49f28..2e26018f 100644 --- a/city_metrix/layers/average_net_building_height.py +++ b/city_metrix/layers/average_net_building_height.py @@ -5,6 +5,7 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection + class AverageNetBuildingHeight(Layer): """ Attributes: @@ -23,8 +24,13 @@ 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, self.spatial_resolution, "average net building height") - .b1) - + + anbh_ic = ee.ImageCollection(anbh) + data = get_image_collection( + anbh_ic, + 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 ab080f51..b04399d9 100644 --- a/city_metrix/layers/built_up_height.py +++ b/city_metrix/layers/built_up_height.py @@ -22,8 +22,15 @@ def get_data(self, bbox): # ANBH is the average height of the built surfaces, USE THIS # AGBH is the amount of built cubic meters per surface unit in the cell # 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.spatial_resolution, "built up height") - return data.built_height + built_height_ic = ee.ImageCollection(built_height) + data = get_image_collection( + built_height_ic, + bbox, + self.spatial_resolution, + "built up height" + ).built_height + + return data diff --git a/city_metrix/layers/esa_world_cover.py b/city_metrix/layers/esa_world_cover.py index a4b0f65c..7e20c31f 100644 --- a/city_metrix/layers/esa_world_cover.py +++ b/city_metrix/layers/esa_world_cover.py @@ -39,19 +39,18 @@ def __init__(self, land_cover_class=None, year=2020, spatial_resolution=10, **kw def get_data(self, bbox): if self.year == 2020: - data = get_image_collection( - ee.ImageCollection("ESA/WorldCover/v100"), - bbox, - self.spatial_resolution, - "ESA world cover" - ).Map + esa_data_ic = ee.ImageCollection("ESA/WorldCover/v100") elif self.year == 2021: - data = get_image_collection( - ee.ImageCollection("ESA/WorldCover/v200"), - bbox, - self.spatial_resolution, - "ESA world cover" - ).Map + esa_data_ic = ee.ImageCollection("ESA/WorldCover/v200") + else: + raise ValueError(f'Specified year ({self.year}) is not currently supported') + + data = get_image_collection( + esa_data_ic, + bbox, + self.spatial_resolution, + "ESA world cover" + ).Map if self.land_cover_class: data = data.where(data == self.land_cover_class.value) diff --git a/city_metrix/layers/glad_lulc.py b/city_metrix/layers/glad_lulc.py index 16d32346..da29076a 100644 --- a/city_metrix/layers/glad_lulc.py +++ b/city_metrix/layers/glad_lulc.py @@ -17,8 +17,9 @@ def __init__(self, year=2020, spatial_resolution=30, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): + lcluc_ic = ee.ImageCollection(ee.Image(f'projects/glad/GLCLU2020/LCLUC_{self.year}')) data = get_image_collection( - ee.ImageCollection(ee.Image(f'projects/glad/GLCLU2020/LCLUC_{self.year}')), + lcluc_ic, bbox, self.spatial_resolution, "GLAD Land Cover" @@ -80,7 +81,8 @@ def __init__(self, year=2020, spatial_resolution=30, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): - simplified_glad = LandCoverSimplifiedGlad(year=self.year, spatial_resolution=self.spatial_resolution).get_data(bbox) + simplified_glad = (LandCoverSimplifiedGlad(year=self.year, spatial_resolution=self.spatial_resolution) + .get_data(bbox)) # Copy the original data data = simplified_glad.copy(deep=True) @@ -108,8 +110,10 @@ def __init__(self, start_year=2000, end_year=2020, spatial_resolution=30, **kwar self.spatial_resolution = spatial_resolution def get_data(self, bbox): - habitat_glad_start = LandCoverHabitatGlad(year=self.start_year, spatial_resolution=self.spatial_resolution).get_data(bbox) - habitat_glad_end = LandCoverHabitatGlad(year=self.end_year, spatial_resolution=self.spatial_resolution).get_data(bbox) + habitat_glad_start = (LandCoverHabitatGlad(year=self.start_year, spatial_resolution=self.spatial_resolution) + .get_data(bbox)) + habitat_glad_end = (LandCoverHabitatGlad(year=self.end_year, spatial_resolution=self.spatial_resolution) + .get_data(bbox)) # Class 01: Became habitat between start year and end year class_01 = ((habitat_glad_start == 0) & (habitat_glad_end == 1)) diff --git a/city_metrix/layers/high_land_surface_temperature.py b/city_metrix/layers/high_land_surface_temperature.py index 5faa2aed..670db62b 100644 --- a/city_metrix/layers/high_land_surface_temperature.py +++ b/city_metrix/layers/high_land_surface_temperature.py @@ -1,12 +1,10 @@ from .landsat_collection_2 import LandsatCollection2 from .land_surface_temperature import LandSurfaceTemperature from .layer import Layer - from shapely.geometry import box import datetime import ee - class HighLandSurfaceTemperature(Layer): """ Attributes: @@ -28,6 +26,7 @@ def get_data(self, bbox): end_date = (hottest_date + datetime.timedelta(days=45)).strftime("%Y-%m-%d") 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 @@ -36,12 +35,17 @@ def get_hottest_date(self, bbox): centroid = box(*bbox).centroid dataset = ee.ImageCollection("ECMWF/ERA5/DAILY") - AirTemperature = (dataset - .filter(ee.Filter.And( - ee.Filter.date(self.start_date, self.end_date), - ee.Filter.bounds(ee.Geometry.BBox(*bbox)))) - .select(['maximum_2m_air_temperature'], ['tasmax']) - ) + + AirTemperature = ( + dataset + .filter( + ee.Filter + .And(ee.Filter.date(self.start_date, self.end_date), + ee.Filter.bounds(ee.Geometry.BBox(*bbox)) + ) + ) + .select(['maximum_2m_air_temperature'], ['tasmax']) + ) # add date as a band to image collection def addDate(image): @@ -56,8 +60,17 @@ def addDate(image): # reduce composite to get the hottest date for centroid of ROI resolution = dataset.first().projection().nominalScale() - hottest_date = str( - ee.Number(hottest.reduceRegion(ee.Reducer.firstNonNull(), ee.Geometry.Point([centroid.x, centroid.y]), resolution).get('date')).getInfo()) + hottest_date = ( + ee.Number( + hottest.reduceRegion(ee.Reducer.firstNonNull(), + ee.Geometry.Point([centroid.x, centroid.y]), + resolution + ).get('date') + ) + .getInfo() + ) # convert to date object - return datetime.datetime.strptime(hottest_date, "%Y%m%d").date() + formated_hottest_data = datetime.datetime.strptime(str(hottest_date), "%Y%m%d").date() + + return formated_hottest_data diff --git a/city_metrix/layers/impervious_surface.py b/city_metrix/layers/impervious_surface.py index ea6772ec..da451598 100644 --- a/city_metrix/layers/impervious_surface.py +++ b/city_metrix/layers/impervious_surface.py @@ -18,12 +18,20 @@ def __init__(self, spatial_resolution=100, **kwargs): def get_data(self, bbox): # load impervious_surface - dataset = ee.ImageCollection(ee.Image("Tsinghua/FROM-GLC/GAIA/v10").gt(0)) # change_year_index is zero if permeable as of 2018 - imperv_surf = ee.ImageCollection(dataset - .filterBounds(ee.Geometry.BBox(*bbox)) - .select('change_year_index') - .sum() - ) - - data = get_image_collection(imperv_surf, bbox, self.spatial_resolution, "imperv surf") - return data.change_year_index + # change_year_index is zero if permeable as of 2018 + impervious_surface = ee.ImageCollection(ee.Image("Tsinghua/FROM-GLC/GAIA/v10").gt(0)) + + imperv_surf_ic = ee.ImageCollection(impervious_surface + .filterBounds(ee.Geometry.BBox(*bbox)) + .select('change_year_index') + .sum() + ) + + data = get_image_collection( + imperv_surf_ic, + bbox, + self.spatial_resolution, + "imperv surf" + ).change_year_index + + return data diff --git a/city_metrix/layers/land_surface_temperature.py b/city_metrix/layers/land_surface_temperature.py index 931cb2e0..48a66e66 100644 --- a/city_metrix/layers/land_surface_temperature.py +++ b/city_metrix/layers/land_surface_temperature.py @@ -1,6 +1,5 @@ from .landsat_collection_2 import LandsatCollection2 from .layer import Layer, get_utm_zone_epsg, get_image_collection - from dask.diagnostics import ProgressBar import ee import xarray @@ -22,6 +21,7 @@ def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resol def get_data(self, bbox): def cloud_mask(image): qa = image.select('QA_PIXEL') + mask = qa.bitwiseAnd(1 << 3).Or(qa.bitwiseAnd(1 << 4)) return image.updateMask(mask.Not()) @@ -30,13 +30,22 @@ def apply_scale_factors(image): return thermal_band l8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2") - l8_st = l8 \ - .select('ST_B10', 'QA_PIXEL') \ - .filter(ee.Filter.date(self.start_date, self.end_date)) \ - .filterBounds(ee.Geometry.BBox(*bbox)) \ - .map(cloud_mask) \ - .map(apply_scale_factors) \ - .reduce(ee.Reducer.mean()) - - data = get_image_collection(ee.ImageCollection(l8_st), bbox, self.spatial_resolution, "LST").ST_B10_mean + + l8_st = (l8 + .select('ST_B10', 'QA_PIXEL') + .filter(ee.Filter.date(self.start_date, self.end_date)) + .filterBounds(ee.Geometry.BBox(*bbox)) + .map(cloud_mask) + .map(apply_scale_factors) + .reduce(ee.Reducer.mean()) + ) + + l8_st_ic = ee.ImageCollection(l8_st) + data = get_image_collection( + l8_st_ic, + bbox, + self.spatial_resolution, + "LST" + ).ST_B10_mean + return data diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 83fa9eb2..635ec257 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -2,6 +2,7 @@ from abc import abstractmethod from typing import Union, Tuple from uuid import uuid4 +# This osgeo import is essential for proper functioning. Do not remove. from osgeo import gdal import ee @@ -324,6 +325,29 @@ def get_stats_funcs(stats_func): return [stats_func] +def set_resampling_for_continuous_raster(image: ee.Image, resampling_method: str, resolution: int, + bbox: tuple[float, float, float, float]): + """ + Function sets the resampling method on the GEE query dictionary for use on continuous raster layers. + GEE only supports bilinear and bicubic interpolation methods. + """ + valid_raster_resampling_methods = ['bilinear', 'bicubic', None] + + if resampling_method not in valid_raster_resampling_methods: + raise ValueError(f"Invalid resampling method ('{resampling_method}'). " + f"Valid methods: {valid_raster_resampling_methods}") + + if resampling_method is None: + data = image + else: + crs = get_utm_zone_epsg(bbox) + data = (image + .resample(resampling_method) + .reproject(crs=crs, scale=resolution)) + + return data + + def get_image_collection( image_collection: ImageCollection, bbox: Tuple[float], @@ -394,4 +418,4 @@ def offset_meters_to_geographic_degrees(decimal_latitude, length_m): lon_degree_offset = abs((length_m / (earth_radius_m * math.cos(math.pi*decimal_latitude/180))) * rad) lat_degree_offset = abs((length_m / earth_radius_m) * rad) - return lon_degree_offset, lat_degree_offset + return lon_degree_offset, lat_degree_offset \ No newline at end of file diff --git a/city_metrix/layers/nasa_dem.py b/city_metrix/layers/nasa_dem.py index b5ac45d8..8f065ccd 100644 --- a/city_metrix/layers/nasa_dem.py +++ b/city_metrix/layers/nasa_dem.py @@ -2,26 +2,43 @@ import xee import xarray as xr -from .layer import Layer, get_image_collection +from .layer import Layer, get_image_collection, set_resampling_for_continuous_raster class NasaDEM(Layer): """ Attributes: spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) + resampling_method: interpolation method used by Google Earth Engine. Default is 'bilinear'. All options are: ('bilinear', 'bicubic', None). """ - def __init__(self, spatial_resolution=30, **kwargs): + def __init__(self, spatial_resolution:int=30, resampling_method:str='bilinear', **kwargs): super().__init__(**kwargs) self.spatial_resolution = spatial_resolution + self.resampling_method = resampling_method + + def get_data(self, bbox: tuple[float, float, float, float]): + nasa_dem = ee.Image("NASA/NASADEM_HGT/001") + + nasa_dem_elev = (ee.ImageCollection(nasa_dem) + .filterBounds(ee.Geometry.BBox(*bbox)) + .select('elevation') + .map(lambda x: + set_resampling_for_continuous_raster(x, + self.resampling_method, + self.spatial_resolution, + bbox + ) + ) + .mean() + ) + + nasa_dem_elev_ic = ee.ImageCollection(nasa_dem_elev) + data = get_image_collection( + nasa_dem_elev_ic, + bbox, + self.spatial_resolution, + "NASA DEM" + ).elevation - def get_data(self, bbox): - dataset = ee.Image("NASA/NASADEM_HGT/001") - nasa_dem = ee.ImageCollection(ee.ImageCollection(dataset) - .filterBounds(ee.Geometry.BBox(*bbox)) - .select('elevation') - .mean() - ) - 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 5efe4a8e..092a30c7 100644 --- a/city_metrix/layers/natural_areas.py +++ b/city_metrix/layers/natural_areas.py @@ -32,7 +32,11 @@ def get_data(self, bbox): } # Perform the reclassification - reclassified_data = reclassify(esa_world_cover, bins=list(reclass_map.keys()), new_values=list(reclass_map.values())) + reclassified_data = reclassify( + esa_world_cover, + bins=list(reclass_map.keys()), + new_values=list(reclass_map.values()) + ) # Apply the original CRS (Coordinate Reference System) reclassified_data = reclassified_data.rio.write_crs(esa_world_cover.rio.crs, inplace=True) diff --git a/city_metrix/layers/ndvi_sentinel2_gee.py b/city_metrix/layers/ndvi_sentinel2_gee.py index fce30219..b31021dc 100644 --- a/city_metrix/layers/ndvi_sentinel2_gee.py +++ b/city_metrix/layers/ndvi_sentinel2_gee.py @@ -1,4 +1,5 @@ import ee + from .layer import Layer, get_image_collection class NdviSentinel2(Layer): @@ -32,6 +33,7 @@ def calculate_ndvi(image): return image.addBands(ndvi) s2 = ee.ImageCollection("COPERNICUS/S2_HARMONIZED") + ndvi = (s2 .filterBounds(ee.Geometry.BBox(*bbox)) .filterDate(start_date, end_date) @@ -41,8 +43,11 @@ 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) + ndvi_mosaic_ic = ee.ImageCollection(ndvi_mosaic) + ndvi_data = get_image_collection( + ndvi_mosaic_ic, + bbox, + self.spatial_resolution, "NDVI" + ).NDVI return ndvi_data diff --git a/city_metrix/layers/smart_surface_lulc.py b/city_metrix/layers/smart_surface_lulc.py index e4759a69..c8e0582e 100644 --- a/city_metrix/layers/smart_surface_lulc.py +++ b/city_metrix/layers/smart_surface_lulc.py @@ -20,15 +20,16 @@ class SmartSurfaceLULC(Layer): - def __init__(self, land_cover_class=None, **kwargs): + def __init__(self, land_cover_class=None, spatial_resolution=10, **kwargs): super().__init__(**kwargs) self.land_cover_class = land_cover_class + self.spatial_resolution = spatial_resolution def get_data(self, bbox): crs = get_utm_zone_epsg(bbox) # ESA world cover - esa_world_cover = EsaWorldCover(year=2021).get_data(bbox) + esa_world_cover = EsaWorldCover(year=2021, spatial_resolution = self.spatial_resolution).get_data(bbox) # ESA reclass and upsample def get_data_esa_reclass(esa_world_cover): reclass_map = { @@ -46,7 +47,11 @@ def get_data_esa_reclass(esa_world_cover): } # Perform the reclassification - reclassified_esa = reclassify(esa_world_cover, bins=list(reclass_map.keys()), new_values=list(reclass_map.values())) + reclassified_esa = reclassify( + esa_world_cover, + bins=list(reclass_map.keys()), + new_values=list(reclass_map.values()) + ) esa_1m = reclassified_esa.rio.reproject( dst_crs=crs, @@ -74,7 +79,8 @@ def get_data_esa_reclass(esa_world_cover): if len(roads_osm) > 0: roads_osm['lanes'] = pd.to_numeric(roads_osm['lanes'], errors='coerce') # Get the average number of lanes per highway class - lanes = (roads_osm.drop(columns='geometry') + lanes = (roads_osm + .drop(columns='geometry') .groupby('highway') # Calculate average and round up .agg(avg_lanes=('lanes', lambda x: np.ceil(np.nanmean(x)) if not np.isnan(x).all() else np.NaN)) @@ -93,11 +99,16 @@ def get_data_esa_reclass(esa_world_cover): # https://nacto.org/publication/urban-street-design-guide/street-design-elements/lane-width/#:~:text=wider%20lane%20widths.-,Lane%20widths%20of%2010%20feet%20are%20appropriate%20in%20urban%20areas,be%20used%20in%20each%20direction # cap is flat to the terminus of the road # join style is mitred so intersections are squared - roads_osm['geometry'] = roads_osm.apply(lambda row: row['geometry'].buffer( - row['lanes'] * 3.048 / 2, - cap_style=CAP_STYLE.flat, - join_style=JOIN_STYLE.mitre), - axis=1 + roads_osm['geometry'] = ( + roads_osm + .apply( + lambda row: row['geometry'] + .buffer( + row['lanes'] * 3.048 / 2, + cap_style=CAP_STYLE.flat, + join_style=JOIN_STYLE.mitre + ),axis=1 + ) ) else: # Add value field (30) @@ -179,14 +190,20 @@ def classify_building(building): if len(buildings) > 0: buildings['Value'] = buildings.apply(classify_building, axis=1) - # Parking parking_osm = OpenStreetMap(osm_class=OpenStreetMapClass.PARKING).get_data(bbox).reset_index() parking_osm['Value'] = 50 - # combine features: open space, water, road, building, parking - feature_df = pd.concat([open_space_osm[['geometry', 'Value']], water_osm[['geometry', 'Value']], roads_osm[['geometry', 'Value']], buildings[['geometry', 'Value']], parking_osm[['geometry', 'Value']]], axis=0) + feature_df = pd.concat( + [open_space_osm[['geometry', 'Value']], + water_osm[['geometry', 'Value']], + roads_osm[['geometry', 'Value']], + buildings[['geometry', 'Value']], + parking_osm[['geometry', 'Value']] + ], axis=0 + ) + # rasterize if feature_df.empty: feature_1m = xr.zeros_like(esa_1m) diff --git a/city_metrix/layers/tree_canopy_height.py b/city_metrix/layers/tree_canopy_height.py index fee1697b..9dbeb524 100644 --- a/city_metrix/layers/tree_canopy_height.py +++ b/city_metrix/layers/tree_canopy_height.py @@ -20,10 +20,19 @@ def __init__(self, spatial_resolution=1, **kwargs): 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") + canopy_ht_img = (canopy_ht + .reduce(ee.Reducer.mean()) + .rename("cover_code") + ) - data = get_image_collection(ee.ImageCollection(canopy_ht), bbox, - self.spatial_resolution, "tree canopy height") + canopy_ht_ic = ee.ImageCollection(canopy_ht_img) + data = get_image_collection( + canopy_ht_ic, + bbox, + self.spatial_resolution, + "tree canopy height" + ).cover_code - return data.cover_code + return data diff --git a/city_metrix/layers/tree_cover.py b/city_metrix/layers/tree_cover.py index 98bc481d..4c1a16e2 100644 --- a/city_metrix/layers/tree_cover.py +++ b/city_metrix/layers/tree_cover.py @@ -5,6 +5,7 @@ import xee import ee + class TreeCover(Layer): """ Merged tropical and nontropical tree cover from WRI @@ -24,11 +25,21 @@ def __init__(self, min_tree_cover=None, max_tree_cover=None, spatial_resolution= def get_data(self, bbox): tropics = ee.ImageCollection('projects/wri-datalab/TropicalTreeCover') - nontropics = ee.ImageCollection('projects/wri-datalab/TTC-nontropics') - 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.spatial_resolution, "tree cover").ttc + non_tropics = ee.ImageCollection('projects/wri-datalab/TTC-nontropics') + + merged_ttc = tropics.merge(non_tropics) + ttc_image = (merged_ttc + .reduce(ee.Reducer.mean()) + .rename('ttc') + ) + + ttc_ic = ee.ImageCollection(ttc_image) + data = get_image_collection( + ttc_ic, + bbox, + self.spatial_resolution, + "tree cover" + ).ttc if self.min_tree_cover is not None: data = data.where(data >= self.min_tree_cover) @@ -36,4 +47,3 @@ def get_data(self, bbox): data = data.where(data <= self.max_tree_cover) return data - diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py index fe69c758..1c81d391 100644 --- a/city_metrix/layers/urban_land_use.py +++ b/city_metrix/layers/urban_land_use.py @@ -19,21 +19,28 @@ def __init__(self, band='lulc', spatial_resolution=5, **kwargs): self.spatial_resolution = spatial_resolution def get_data(self, bbox): - dataset = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1") + ulu = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1") + # ImageCollection didn't cover the globe - if dataset.filterBounds(ee.Geometry.BBox(*bbox)).size().getInfo() == 0: - ulu = ee.ImageCollection(ee.Image.constant(0) + if ulu.filterBounds(ee.Geometry.BBox(*bbox)).size().getInfo() == 0: + ulu_ic = ee.ImageCollection(ee.Image + .constant(0) .clip(ee.Geometry.BBox(*bbox)) .rename('lulc') ) else: - ulu = ee.ImageCollection(dataset + ulu_ic = ee.ImageCollection(ulu .filterBounds(ee.Geometry.BBox(*bbox)) .select(self.band) .reduce(ee.Reducer.firstNonNull()) .rename('lulc') ) - data = get_image_collection(ulu, bbox, self.spatial_resolution, "urban land use").lulc + data = get_image_collection( + ulu_ic, + 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 30fb6d8d..88a5ad09 100644 --- a/city_metrix/layers/world_pop.py +++ b/city_metrix/layers/world_pop.py @@ -27,23 +27,44 @@ def get_data(self, bbox): if not self.agesex_classes: # total population dataset = ee.ImageCollection('WorldPop/GP/100m/pop') - world_pop = ee.ImageCollection(dataset - .filterBounds(ee.Geometry.BBox(*bbox)) - .filter(ee.Filter.inList('year', [self.year])) - .select('population') - .mean() - ) - data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop").population + world_pop_ic = ee.ImageCollection( + dataset + .filterBounds(ee.Geometry.BBox(*bbox)) + .filter(ee.Filter.inList('year', [self.year])) + .select('population') + .mean() + ) + + data = get_image_collection( + world_pop_ic, + bbox, + self.spatial_resolution, + "world pop" + ).population else: # sum population for selected age-sex groups - dataset = ee.ImageCollection('WorldPop/GP/100m/pop_age_sex') - world_pop = dataset.filterBounds(ee.Geometry.BBox(*bbox))\ - .filter(ee.Filter.inList('year', [self.year]))\ - .select(self.agesex_classes)\ - .mean() - world_pop = ee.ImageCollection(world_pop.reduce(ee.Reducer.sum()).rename('sum_age_sex_group')) - data = get_image_collection(world_pop, bbox, self.spatial_resolution, "world pop age sex").sum_age_sex_group + world_pop_age_sex = ee.ImageCollection('WorldPop/GP/100m/pop_age_sex') + + world_pop_age_sex_year = (world_pop_age_sex + .filterBounds(ee.Geometry.BBox(*bbox)) + .filter(ee.Filter.inList('year', [self.year])) + .select(self.agesex_classes) + .mean() + ) + + world_pop_group_ic = ee.ImageCollection( + world_pop_age_sex_year + .reduce(ee.Reducer.sum()) + .rename('sum_age_sex_group') + ) + + data = get_image_collection( + world_pop_group_ic, + bbox, + self.spatial_resolution, + "world pop age sex" + ).sum_age_sex_group return data diff --git a/tests/resources/bbox_constants.py b/tests/resources/bbox_constants.py index 789ab48f..57bd552b 100644 --- a/tests/resources/bbox_constants.py +++ b/tests/resources/bbox_constants.py @@ -24,3 +24,12 @@ -38.39993,-12.93239 ) +BBOX_NLD_AMSTERDAM_TEST = ( + 4.9012,52.3720, + 4.9083,52.3752 +) + +BBOX_NLD_AMSTERDAM_LARGE_TEST = ( + 4.884629880473071,52.34146514406914, + 4.914180290924863,52.359560786247165 +) diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/README.md b/tests/resources/layer_dumps/README.md similarity index 100% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/README.md rename to tests/resources/layer_dumps/README.md diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/__init__.py b/tests/resources/layer_dumps/__init__.py similarity index 100% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/__init__.py rename to tests/resources/layer_dumps/__init__.py diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py b/tests/resources/layer_dumps/conftest.py similarity index 84% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py rename to tests/resources/layer_dumps/conftest.py index d72876d7..9a5dd9a6 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py +++ b/tests/resources/layer_dumps/conftest.py @@ -4,7 +4,8 @@ import shutil from collections import namedtuple -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_NLD_AMSTERDAM_TEST, \ + BBOX_NLD_AMSTERDAM_LARGE_TEST 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 @@ -19,23 +20,20 @@ # 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 +# BBOX = BBOX_NLD_AMSTERDAM_TEST +# BBOX = BBOX_NLD_AMSTERDAM_LARGE_TEST # 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' if RUN_DUMPS is True: source_folder = os.path.dirname(__file__) target_folder = get_target_folder_path() - create_target_folder(target_folder, True) + create_target_folder(target_folder, False) - source_qgis_file = os.path.join(source_folder, qgis_project_file) - target_qgis_file = os.path.join(target_folder, qgis_project_file) - shutil.copyfile(source_qgis_file, target_qgis_file) - - print("\n\033[93m QGIS project file and layer files written to folder %s.\033[0m\n" % target_folder) + print("\n\033[93m Layer files written to folder %s.\033[0m\n" % target_folder) @pytest.fixture def target_folder(): diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py b/tests/resources/layer_dumps/test_write_all_layers.py similarity index 81% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py rename to tests/resources/layer_dumps/test_write_all_layers.py index 9c382116..2be8c4df 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_to_qgis_files.py +++ b/tests/resources/layer_dumps/test_write_all_layers.py @@ -2,27 +2,7 @@ # Execution configuration is specified in the conftest file import pytest -from city_metrix.layers import ( - Albedo, - AlosDSM, - AverageNetBuildingHeight, - EsaWorldCover, - HighLandSurfaceTemperature, - LandsatCollection2, - LandSurfaceTemperature, - NasaDEM, - NaturalAreas, - OpenBuildings, - OpenStreetMap, - OvertureBuildings, - Sentinel2Level2, - NdviSentinel2, - SmartSurfaceLULC, - TreeCanopyHeight, - TreeCover, - UrbanLandUse, - WorldPop, Layer, ImperviousSurface -) +from city_metrix.layers import * from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated, get_file_count_in_folder from ...tools.general_tools import get_class_default_spatial_resolution @@ -34,30 +14,6 @@ def test_write_albedo(target_folder, bbox_info, target_spatial_resolution_multip 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_albedo_tiled_unbuffered(target_folder, bbox_info, target_spatial_resolution_multiplier): - file_path = prep_output_path(target_folder, 'albedo_tiled.tif') - 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=0.01, buffer_size=None)) - file_count = get_file_count_in_folder(file_path) - - expected_file_count = 5 # includes 4 tiles and one geojson file - assert file_count == expected_file_count - -@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') -def test_write_albedo_tiled_buffered(target_folder, bbox_info, target_spatial_resolution_multiplier): - buffer_degrees = 0.001 - file_path = prep_output_path(target_folder, 'albedo_tiled_buffered.tif') - 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=0.01, buffer_size=buffer_degrees)) - file_count = get_file_count_in_folder(file_path) - - expected_file_count = 6 # includes 4 tiles and two geojson files - assert file_count == expected_file_count - - @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') def test_write_alos_dsm(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'alos_dsm.tif') @@ -72,6 +28,13 @@ def test_write_average_net_building_height(target_folder, bbox_info, target_spat 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_built_up_height(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, 'built_up_height.tif') + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(BuiltUpHeight()) + BuiltUpHeight(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, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'esa_world_cover.tif') @@ -79,6 +42,13 @@ def test_write_esa_world_cover(target_folder, bbox_info, target_spatial_resoluti 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_glad_lulc(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, 'glad_lulc.tif') + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(LandCoverGlad()) + LandCoverGlad(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, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'high_land_surface_temperature.tif') @@ -90,14 +60,7 @@ def test_write_high_land_surface_temperature(target_folder, bbox_info, target_sp def test_write_impervious_surface(target_folder, bbox_info, target_spatial_resolution_multiplier): file_path = prep_output_path(target_folder, 'impervious_surface.tif') target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(ImperviousSurface()) - LandSurfaceTemperature(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, target_spatial_resolution_multiplier): - file_path = prep_output_path(target_folder, 'land_surface_temperature.tif') - 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) + ImperviousSurface(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 @@ -158,9 +121,9 @@ def test_write_overture_buildings(target_folder, bbox_info, target_spatial_resol @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') 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) + target_resolution = target_spatial_resolution_multiplier * get_class_default_spatial_resolution(SmartSurfaceLULC()) + SmartSurfaceLULC(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') diff --git a/tests/resources/layer_dumps/test_write_layers_other.py b/tests/resources/layer_dumps/test_write_layers_other.py new file mode 100644 index 00000000..d4b0a23b --- /dev/null +++ b/tests/resources/layer_dumps/test_write_layers_other.py @@ -0,0 +1,31 @@ +# This code is mostly intended for manual execution +# Execution configuration is specified in the conftest file +import pytest + +from city_metrix.layers import * +from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated, get_file_count_in_folder +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_tiled_unbuffered(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, 'albedo_tiled.tif') + 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=0.01, buffer_size=None)) + file_count = get_file_count_in_folder(file_path) + + expected_file_count = 5 # includes 4 tiles and one geojson file + assert file_count == expected_file_count + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_albedo_tiled_buffered(target_folder, bbox_info, target_spatial_resolution_multiplier): + buffer_degrees = 0.001 + file_path = prep_output_path(target_folder, 'albedo_tiled_buffered.tif') + 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=0.01, buffer_size=buffer_degrees)) + file_count = get_file_count_in_folder(file_path) + + expected_file_count = 6 # includes 4 tiles and two geojson files + assert file_count == expected_file_count \ No newline at end of file diff --git a/tests/resources/layer_dumps/test_write_layers_using_fixed_resolution.py b/tests/resources/layer_dumps/test_write_layers_using_fixed_resolution.py new file mode 100644 index 00000000..bf199916 --- /dev/null +++ b/tests/resources/layer_dumps/test_write_layers_using_fixed_resolution.py @@ -0,0 +1,116 @@ +# This code is mostly intended for manual execution +# Execution configuration is specified in the conftest file +import pytest + +from city_metrix.layers import * +from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated + +TARGET_RESOLUTION = 5 +TARGET_RESAMPLING_METHOD = 'bilinear' +# TARGET_RESAMPLING_METHOD = None + +@pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False') +def test_write_albedo_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'albedo_{TARGET_RESOLUTION}m.tif') + (Albedo(spatial_resolution=TARGET_RESOLUTION, resampling_method=TARGET_RESAMPLING_METHOD) + .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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'alos_dsm_{TARGET_RESOLUTION}m.tif') + (AlosDSM(spatial_resolution=TARGET_RESOLUTION, resampling_method=TARGET_RESAMPLING_METHOD) + .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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'average_net_building_height_{TARGET_RESOLUTION}m.tif') + 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_built_up_height_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'built_up_height_{TARGET_RESOLUTION}m.tif') + BuiltUpHeight(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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'esa_world_cover_{TARGET_RESOLUTION}m.tif') + 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_glad_lulc_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'glad_lulc_{TARGET_RESOLUTION}m.tif') + LandCoverGlad(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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'high_land_surface_temperature_{TARGET_RESOLUTION}m.tif') + 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_impervious_surface_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'impervious_surface_{TARGET_RESOLUTION}m.tif') + ImperviousSurface(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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'land_surface_temperature_{TARGET_RESOLUTION}m.tif') + LandSurfaceTemperature(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_nasa_dem_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'nasa_dem_{TARGET_RESOLUTION}m.tif') + (NasaDEM(spatial_resolution=TARGET_RESOLUTION, resampling_method=TARGET_RESAMPLING_METHOD) + .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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'natural_areas_{TARGET_RESOLUTION}m.tif') + 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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'ndvi_sentinel2_gee_{TARGET_RESOLUTION}m.tif') + 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_smart_surface_lulc_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'smart_surface_lulc_{TARGET_RESOLUTION}m.tif') + SmartSurfaceLULC(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_canopy_height_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'tree_canopy_height_{TARGET_RESOLUTION}m.tif') + 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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'tree_cover_{TARGET_RESOLUTION}m.tif') + 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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'urban_land_use_{TARGET_RESOLUTION}m.tif') + 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_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier): + file_path = prep_output_path(target_folder, f'world_pop_{TARGET_RESOLUTION}m.tif') + 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/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz b/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz deleted file mode 100644 index 1cc53ce2..00000000 Binary files a/tests/resources/layer_dumps_for_br_lauro_de_freitas/layers_for_br_lauro_de_freitas.qgz and /dev/null differ diff --git a/tests/test_layer_metrics.py b/tests/test_layer_metrics.py index 0002f580..94d558ae 100644 --- a/tests/test_layer_metrics.py +++ b/tests/test_layer_metrics.py @@ -47,8 +47,23 @@ def test_read_image_collection_scale(): pytest.approx(expected_y_size, rel=EE_IMAGE_DIMENSION_TOLERANCE) == actual_y_size ) -def test_albedo_metrics(): - data = Albedo().get_data(BBOX) +def test_albedo_metrics_default_resampling(): + # Default resampling_method is bilinear + data = Albedo(spatial_resolution=10).get_data(BBOX) + + # Bounding values + expected_min_value = _convert_fraction_to_rounded_percent(0.03) + expected_max_value = _convert_fraction_to_rounded_percent(0.34) + actual_min_value = _convert_fraction_to_rounded_percent(data.values.min()) + actual_max_value = _convert_fraction_to_rounded_percent(data.values.max()) + + # Value range + assert expected_min_value == actual_min_value + assert expected_max_value == actual_max_value + + +def test_albedo_metrics_no_resampling(): + data = Albedo(spatial_resolution=10, resampling_method= None).get_data(BBOX) # Bounding values expected_min_value = _convert_fraction_to_rounded_percent(0.03) @@ -62,7 +77,7 @@ def test_albedo_metrics(): def test_alos_dsm_metrics(): - data = AlosDSM().get_data(BBOX) + data = AlosDSM(resampling_method=None).get_data(BBOX) # Bounding values expected_min_value = 16 diff --git a/tests/test_layer_parameters.py b/tests/test_layer_parameters.py index 3fe415be..359bb287 100644 --- a/tests/test_layer_parameters.py +++ b/tests/test_layer_parameters.py @@ -1,4 +1,5 @@ import pytest +import xarray as xr import numpy as np from skimage.metrics import structural_similarity as ssim from pyproj import CRS @@ -286,7 +287,10 @@ def _get_sample_data(class_instance, bbox, downsize_factor): 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'] + if isinstance(image_data, xr.DataArray): + crs_string = image_data.rio.crs.wkt + else: + crs_string = image_data.crs crs = CRS.from_string(crs_string) return crs