diff --git a/README.md b/README.md index 56ae15e..4af03c1 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 ac6b6b4..9b35945 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -1,6 +1,9 @@ import ee +import xarray +from dask.diagnostics import ProgressBar + +from .layer import Layer, get_utm_zone_epsg, get_image_collection, set_resampling_for_continuous_raster -from .layer import Layer, get_image_collection class Albedo(Layer): """ @@ -8,16 +11,19 @@ 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 """ MAX_CLOUD_PROB = 30 S2_ALBEDO_EQN = '((B*Bw)+(G*Gw)+(R*Rw)+(NIR*NIRw)+(SWIR1*SWIR1w)+(SWIR2*SWIR2w))' - 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 ## METHODS @@ -109,7 +115,17 @@ def calc_s2_albedo(image): # S2 MOSAIC AND ALBEDO dataset = self.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()) + + 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( diff --git a/city_metrix/layers/alos_dsm.py b/city_metrix/layers/alos_dsm.py index c22df82..f9601eb 100644 --- a/city_metrix/layers/alos_dsm.py +++ b/city_metrix/layers/alos_dsm.py @@ -1,26 +1,40 @@ import ee +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): + 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') - .mean() - ) + 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() + ) + data = get_image_collection( alos_dsm_ic, diff --git a/city_metrix/layers/average_net_building_height.py b/city_metrix/layers/average_net_building_height.py index 11799cc..2e26018 100644 --- a/city_metrix/layers/average_net_building_height.py +++ b/city_metrix/layers/average_net_building_height.py @@ -1,6 +1,9 @@ +from dask.diagnostics import ProgressBar +import xarray as xr +import xee import ee -from .layer import Layer, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection class AverageNetBuildingHeight(Layer): diff --git a/city_metrix/layers/built_up_height.py b/city_metrix/layers/built_up_height.py index aef268e..b04399d 100644 --- a/city_metrix/layers/built_up_height.py +++ b/city_metrix/layers/built_up_height.py @@ -1,6 +1,9 @@ +from dask.diagnostics import ProgressBar +import xarray as xr +import xee import ee -from .layer import Layer, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection class BuiltUpHeight(Layer): diff --git a/city_metrix/layers/esa_world_cover.py b/city_metrix/layers/esa_world_cover.py index c214781..7e20c31 100644 --- a/city_metrix/layers/esa_world_cover.py +++ b/city_metrix/layers/esa_world_cover.py @@ -1,7 +1,9 @@ +from dask.diagnostics import ProgressBar +from enum import Enum +import xarray as xr import ee -from enum import Enum -from .layer import Layer, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection class EsaWorldCoverClass(Enum): diff --git a/city_metrix/layers/glad_lulc.py b/city_metrix/layers/glad_lulc.py index 11afe31..da29076 100644 --- a/city_metrix/layers/glad_lulc.py +++ b/city_metrix/layers/glad_lulc.py @@ -1,7 +1,7 @@ import xarray as xr import ee -from .layer import Layer, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection class LandCoverGlad(Layer): diff --git a/city_metrix/layers/high_land_surface_temperature.py b/city_metrix/layers/high_land_surface_temperature.py index 9cbf5ea..670db62 100644 --- a/city_metrix/layers/high_land_surface_temperature.py +++ b/city_metrix/layers/high_land_surface_temperature.py @@ -1,9 +1,9 @@ -import datetime -import ee - -from shapely.geometry import box +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): """ diff --git a/city_metrix/layers/impervious_surface.py b/city_metrix/layers/impervious_surface.py index 9cf017d..da45159 100644 --- a/city_metrix/layers/impervious_surface.py +++ b/city_metrix/layers/impervious_surface.py @@ -1,6 +1,9 @@ +from dask.diagnostics import ProgressBar +import xarray as xr +import xee import ee -from .layer import Layer, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection class ImperviousSurface(Layer): diff --git a/city_metrix/layers/land_surface_temperature.py b/city_metrix/layers/land_surface_temperature.py index 93888ce..48a66e6 100644 --- a/city_metrix/layers/land_surface_temperature.py +++ b/city_metrix/layers/land_surface_temperature.py @@ -1,6 +1,8 @@ +from .landsat_collection_2 import LandsatCollection2 +from .layer import Layer, get_utm_zone_epsg, get_image_collection +from dask.diagnostics import ProgressBar import ee - -from .layer import Layer, get_image_collection +import xarray class LandSurfaceTemperature(Layer): """ diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index 3d2829f..635ec25 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -2,6 +2,8 @@ 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 import boto3 @@ -323,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], @@ -393,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 d3840d3..8f065cc 100644 --- a/city_metrix/layers/nasa_dem.py +++ b/city_metrix/layers/nasa_dem.py @@ -1,24 +1,35 @@ import ee +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): + 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() ) @@ -29,5 +40,5 @@ def get_data(self, 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 fdf499b..092a30c 100644 --- a/city_metrix/layers/natural_areas.py +++ b/city_metrix/layers/natural_areas.py @@ -1,3 +1,4 @@ +import xarray as xr from xrspatial.classify import reclassify from .layer import Layer diff --git a/city_metrix/layers/smart_surface_lulc.py b/city_metrix/layers/smart_surface_lulc.py index 2e0f083..c8e0582 100644 --- a/city_metrix/layers/smart_surface_lulc.py +++ b/city_metrix/layers/smart_surface_lulc.py @@ -1,7 +1,9 @@ import xarray as xr import numpy as np import pandas as pd +import geopandas as gpd from shapely.geometry import CAP_STYLE, JOIN_STYLE +from shapely.geometry import box from exactextract import exact_extract from geocube.api.core import make_geocube import warnings diff --git a/city_metrix/layers/tree_canopy_height.py b/city_metrix/layers/tree_canopy_height.py index f499552..9dbeb52 100644 --- a/city_metrix/layers/tree_canopy_height.py +++ b/city_metrix/layers/tree_canopy_height.py @@ -1,6 +1,9 @@ -import ee +from .layer import Layer, get_utm_zone_epsg, get_image_collection -from .layer import Layer, get_image_collection +from dask.diagnostics import ProgressBar +import xarray as xr +import xee +import ee class TreeCanopyHeight(Layer): """ diff --git a/city_metrix/layers/tree_cover.py b/city_metrix/layers/tree_cover.py index 52fc8d6..4c1a16e 100644 --- a/city_metrix/layers/tree_cover.py +++ b/city_metrix/layers/tree_cover.py @@ -1,6 +1,10 @@ +from .layer import Layer, get_utm_zone_epsg, get_image_collection + +from dask.diagnostics import ProgressBar +import xarray as xr +import xee import ee -from .layer import Layer, get_image_collection class TreeCover(Layer): """ diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py index 02c737c..1c81d39 100644 --- a/city_metrix/layers/urban_land_use.py +++ b/city_metrix/layers/urban_land_use.py @@ -1,6 +1,9 @@ +from dask.diagnostics import ProgressBar +import xarray as xr +import xee import ee -from .layer import Layer, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection class UrbanLandUse(Layer): diff --git a/city_metrix/layers/world_pop.py b/city_metrix/layers/world_pop.py index 700010a..88a5ad0 100644 --- a/city_metrix/layers/world_pop.py +++ b/city_metrix/layers/world_pop.py @@ -1,6 +1,9 @@ +from dask.diagnostics import ProgressBar +import xarray as xr +import xee import ee -from .layer import Layer, get_image_collection +from .layer import Layer, get_utm_zone_epsg, get_image_collection class WorldPop(Layer): diff --git a/tests/resources/bbox_constants.py b/tests/resources/bbox_constants.py index 789ab48..57bd552 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 94% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py rename to tests/resources/layer_dumps/conftest.py index 6b33b51..9a5dd9a 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,6 +20,8 @@ # 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 diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_all_layers.py b/tests/resources/layer_dumps/test_write_all_layers.py similarity index 100% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_all_layers.py rename to tests/resources/layer_dumps/test_write_all_layers.py diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_other.py b/tests/resources/layer_dumps/test_write_layers_other.py similarity index 100% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_other.py rename to tests/resources/layer_dumps/test_write_layers_other.py diff --git a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_using_fixed_resolution.py b/tests/resources/layer_dumps/test_write_layers_using_fixed_resolution.py similarity index 93% rename from tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_using_fixed_resolution.py rename to tests/resources/layer_dumps/test_write_layers_using_fixed_resolution.py index 8964b86..bf19991 100644 --- a/tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_using_fixed_resolution.py +++ b/tests/resources/layer_dumps/test_write_layers_using_fixed_resolution.py @@ -3,20 +3,24 @@ 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 .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).write(bbox_info.bounds, file_path, tile_degrees=None) + (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).write(bbox_info.bounds, file_path, tile_degrees=None) + (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') @@ -64,7 +68,8 @@ def test_write_land_surface_temperature_fixed_res(target_folder, bbox_info, targ @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).write(bbox_info.bounds, file_path, tile_degrees=None) + (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') diff --git a/tests/test_layer_metrics.py b/tests/test_layer_metrics.py index 0002f58..94d558a 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