From 6a058441784586544a1073d778c5553a066f43d5 Mon Sep 17 00:00:00 2001 From: weiqi-tori Date: Wed, 23 Oct 2024 19:07:18 +0800 Subject: [PATCH 1/6] adjust albedo layer and add veg_water_map layer --- city_metrix/layers/albedo.py | 91 ++++----- city_metrix/layers/vegetation_water_map.py | 223 +++++++++++++++++++++ 2 files changed, 257 insertions(+), 57 deletions(-) create mode 100644 city_metrix/layers/vegetation_water_map.py diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py index 341786b2..38983397 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -4,6 +4,7 @@ from .layer import Layer, get_utm_zone_epsg, get_image_collection + class Albedo(Layer): """ Attributes: @@ -12,6 +13,11 @@ class Albedo(Layer): spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) threshold: threshold value for filtering the retrieval """ + S2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") + S2C = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY") + + 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): super().__init__(**kwargs) @@ -20,56 +26,31 @@ def __init__(self, start_date="2021-01-01", end_date="2022-01-01", spatial_resol self.spatial_resolution = spatial_resolution self.threshold = threshold - def get_data(self, bbox): - S2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") - S2C = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY") - - MAX_CLOUD_PROB = 30 - S2_ALBEDO_EQN = '((B*Bw)+(G*Gw)+(R*Rw)+(NIR*NIRw)+(SWIR1*SWIR1w)+(SWIR2*SWIR2w))' - - ## METHODS - - ## get cloudmasked image collection - - 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') - 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) - - 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) - - def get_masked_s2_collection(roi, start, end): - criteria = (ee.Filter.And( - ee.Filter.date(start, end), - ee.Filter.bounds(roi) - )) - 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'}) - })) - - 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) + # get cloudmasked image collection + def mask_clouds_and_rescale(self, im): + clouds = ee.Image(im.get('cloud_mask')).select('probability') + return im.updateMask(clouds.lt(self.MAX_CLOUD_PROB)).divide(10000) + + def get_masked_s2_collection(self, roi, start, end): + criteria = (ee.Filter.And( + ee.Filter.date(start, end), + ee.Filter.bounds(roi) + )) + # .select('B2','B3','B4','B8','B11','B12') + s2 = self.S2.filter(criteria) + s2c = self.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=s2_with_clouds.limit(image_limit,'nb_cloudy_pixels') + s2_with_clouds = ee.ImageCollection(s2_with_clouds).map( + self.mask_clouds_and_rescale) # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE') + return ee.ImageCollection(s2_with_clouds) + def get_data(self, bbox): # calculate albedo for images # weights derived from @@ -89,20 +70,16 @@ def calc_s2_albedo(image): 'SWIR1': image.select('B11'), 'SWIR2': image.select('B12') } - return image.expression(S2_ALBEDO_EQN, config).double().rename('albedo') + return image.expression(self.S2_ALBEDO_EQN, config).double().rename('albedo') - ## S2 MOSAIC AND ALBEDO - dataset = get_masked_s2_collection(ee.Geometry.BBox(*bbox), self.start_date, self.end_date) + # 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()) - data = (get_image_collection( - ee.ImageCollection(albedo_mean), bbox, self.spatial_resolution, "albedo") - .albedo_mean) + data = (get_image_collection(ee.ImageCollection(albedo_mean), 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/vegetation_water_map.py b/city_metrix/layers/vegetation_water_map.py new file mode 100644 index 00000000..b936de4e --- /dev/null +++ b/city_metrix/layers/vegetation_water_map.py @@ -0,0 +1,223 @@ +import xarray as xr +import ee + +from .layer import Layer, get_image_collection +from .albedo import Albedo + + +class VegetationWaterMap(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) + """ + + def __init__(self, start_date="2018-01-01", end_date="2022-12-31", 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): + NDVIthreshold = 0.4 # decimal + NDWIthreshold = 0.3 # decimal + # half the value of the p-value threshold to be used in the significance test. + halfpvalue = 0.025 + + # annual image collections and images + def AnnualIC(ic, year): + def addNDVI(image): + ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI') + return image.addBands(ndvi) + def addNDWI(image): + ndwi = image.normalizedDifference(['B3', 'B8']).rename('NDWI') + return image.addBands(ndwi) + def addyear(image): + return image.set("year", year) + ImgColl = ic.filter(ee.Filter.stringStartsWith('system:index', year)).map(addNDVI).map(addNDWI).map(addyear) + return ImgColl + + def AnnualImgGreen(ic, year): + greenest = ic.select('NDVI').reduce(ee.Reducer.median()).rename('NDVI').addBands(ee.Image(year).rename('time_start')).float() + bluest = ic.qualityMosaic('NDWI').select('NDWI').addBands(ee.Image(year).rename('time_start')).float() + bluestThres = bluest.updateMask(bluest.select('NDWI').gte(NDWIthreshold)) + greenestnowater = greenest.updateMask(bluestThres.select('NDWI').unmask().Not()) + return greenestnowater + def AnnualImgWater(ic, year): + bluest = ic.select('NDWI').reduce(ee.Reducer.median()).rename('NDWI').addBands(ee.Image(year).rename('time_start')).float() + return bluest + def AnnualImgWatermask(ic, year): + bluest = ic.qualityMosaic('NDWI').select('NDWI').addBands(ee.Image(year).rename('time_start')).float() + bluestThres = bluest.updateMask(bluest.select('NDWI').gte(NDWIthreshold)) + return bluestThres + def AnnualImgGreenmask(ic, year): + greenest = ic.qualityMosaic('NDVI').select('NDVI').addBands(ee.Image(year).rename('time_start')).float() + greenestThres = greenest.updateMask(greenest.select('NDVI').gte(NDVIthreshold)) + return greenestThres + + # Functions for significance test + # https://developers.google.com/earth-engine/tutorials/community/nonparametric-trends + # https://code.earthengine.google.com/bce3dc00c56df4246c5d32f5fcccf5c7 + # //////////////////////////////////////singificance test/////////////////////////////////////////////////////////////// + def senum(lfx): + def senumwrap(Im): + esty = lfx.select('scale').multiply(Im.select('time_start')).add(lfx.select('offset')) + diff = Im.select('NDVI').subtract(esty) + pow = diff.multiply(diff) + return (pow) + return (senumwrap) + + def sedenom(mosaicmeanx): + def sedenomwrap(Im): + diff = Im.select('time_start').subtract(mosaicmeanx.select('time_start')) + pow = diff.multiply(diff) + return (pow) + return (sedenomwrap) + + # https://en.wikipedia.org/wiki/Error_function#Cumulative_distribution_function + def eeCdf(t): + return ee.Image(0.5).multiply(ee.Image(1).add(ee.Image(t).divide(ee.Image(2).sqrt()).erf())) + + def invCdf(p): + return ee.Image(2).sqrt().multiply(ee.Image(p).multiply(2).subtract(1).erfInv()) + + # /////green normalized Difference/////// + def significanceGreen(gtrend, glinearfit, SampleNumber): + mosaicmean = gtrend.mean() + mosaicNum = gtrend.map(senum(glinearfit)).sum() + mosaicDenom = gtrend.map(sedenom(mosaicmean)).sum() + StdDev = mosaicNum.divide(mosaicDenom) + # Standard Error - from SampleNumber samples + se = StdDev.divide(SampleNumber).sqrt() + # Test Statistic + gT = glinearfit.select('scale').divide(se) + # Compute P-values. + gP = ee.Image(1).subtract(eeCdf(gT.abs())) + return gP + + # /////water normalized Difference/////// + def significanceWater(wtrend, wlinearfit, SampleNumber): + def senum(lfx): + def senumwrap(Im): + esty = lfx.select('scale').multiply(Im.select('time_start')).add(lfx.select('offset')) + diff = Im.select('NDWI').subtract(esty) + pow = diff.multiply(diff) + return (pow) + return (senumwrap) + + mosaicmean = wtrend.mean() + mosaicNum = wtrend.map(senum(wlinearfit)).sum() + mosaicDenom = wtrend.map(sedenom(mosaicmean)).sum() + StdDev = mosaicNum.divide(mosaicDenom) + # Standard Error - from SampleNumber samples + se = StdDev.divide(SampleNumber).sqrt() + # Test Statistic + wT = wlinearfit.select('scale').divide(se) + # Compute P-values. + wP = ee.Image(1).subtract(eeCdf(wT.abs())) + return wP + + # function to generate vegetation and water trend and change maps + + def get_map_vegwaterchange(IC): + gwic2019 = AnnualIC(IC, '2019') + gwic2020 = AnnualIC(IC, '2020') + gwic2021 = AnnualIC(IC, '2021') + gwic2022 = AnnualIC(IC, '2022') + gwic = ((gwic2019).merge(gwic2020).merge(gwic2021).merge(gwic2022)) + + d = {} + for i in range(2019, 2023): + # https://stackoverflow.com/questions/6181935/how-do-you-create-different-variable-names-while-in-a-loop + filteredgwic = gwic.filter(ee.Filter.eq("year", str(i))) + d["g{0}".format(i)] = AnnualImgGreen(filteredgwic, i) + d["w{0}".format(i)] = AnnualImgWater(filteredgwic, i) + d["w{0}mask".format(i)] = AnnualImgWatermask(filteredgwic, i) + d["g{0}mask".format(i)] = AnnualImgGreenmask(filteredgwic, i) + + gtrend = ee.ImageCollection.fromImages([d["g2019"], d["g2020"], d["g2021"], d["g2022"]]) + wtrend = ee.ImageCollection.fromImages([d["w2019"], d["w2020"], d["w2021"], d["w2022"]]) + wanyyear = (ee.Image( + d["w2019mask"].select('NDWI')) + .blend(d["w2020mask"].select('NDWI')) + .blend(d["w2021mask"].select('NDWI')) + .blend(d["w2022mask"].select('NDWI')) + ) + ganyyear = (ee.Image( + d["g2019mask"].select('NDVI')) + .blend(d["g2020mask"].select('NDVI')) + .blend(d["g2021mask"].select('NDVI')) + .blend(d["g2022mask"].select('NDVI')) + ) + startyear = ee.Number(2019) + endyear = ee.Number(2022) + sampleNumber = (endyear.subtract(startyear)).add(1) + gstartmask = AnnualImgGreenmask(gwic.filter(ee.Filter.eq("year", '2019')), startyear) + wstartmask = AnnualImgWatermask(gwic.filter(ee.Filter.eq("year", '2019')), startyear) + gendmask = AnnualImgGreenmask(gwic.filter(ee.Filter.eq("year", '2022')), endyear) + wendmask = AnnualImgWatermask(gwic.filter(ee.Filter.eq("year", '2022')), endyear) + + # 0.05 # minimum threshold scale value from linear fit trend line to be considered a vegetation change. + minSlopeVeg = ee.Number(0.1) + # 0.05 # minimum threshold scale value from linear fit trend line to be considered a water change. + minSlopeWater = ee.Number(0.1) + + # Create linear fit trend line for years of NDVI data + # https://developers.google.com/earth-engine/guides/reducers_regression#linearfit + glinearfit = ee.Image(gtrend.select(['time_start', 'NDVI']).reduce(ee.Reducer.linearFit())) + + # apply significance test + # Pixels that can have the null hypothesis (there is no trend) rejected. + # Specifically, if the true trend is zero, there would be less than 5% (double "halfpvalue") + # chance of randomly obtaining the observed result (that there is a trend). + gsignifMask = significanceGreen(gtrend, glinearfit, sampleNumber).lte(halfpvalue) + glinearfit = glinearfit.updateMask(gsignifMask) + + # Annual slope value for each pixel above threshold. If interested in value for timeperiod, can use .multiply(ee.Image(6)). Can also mask based on offset to limit based on starting vegetation level. + glfLimit = (glinearfit.select('scale') # .multiply(ee.Image(6)) + .updateMask(glinearfit.select('scale').gte(minSlopeVeg).Or(glinearfit.select('scale').lte(ee.Number(0).subtract(minSlopeVeg)))) + ) + greenanyyearMask = ee.Image(0).where(ganyyear.neq(0), 1) + + glfLimitanyyeargreen = glfLimit.updateMask(greenanyyearMask) + glfLimitanyyeargreenLoss = glfLimitanyyeargreen.updateMask(glfLimitanyyeargreen.lt(0)) + glfLimitanyyeargreenGain = glfLimitanyyeargreen.updateMask(glfLimitanyyeargreen.gt(0)) + + # Create linear fit trend line for years of NDWI data + wlinearfit = ee.Image(wtrend.select(['time_start', 'NDWI']).reduce(ee.Reducer.linearFit())) + + # apply significance test + # Pixels that can have the null hypothesis (there is no trend) rejected. + # Specifically, if the true trend is zero, there would be less than 5% (double "halfpvalue") + # chance of randomly obtaining the observed result (that there is a trend). + wsignifMask = significanceWater(wtrend, wlinearfit, sampleNumber).lte(halfpvalue) + wlinearfit = wlinearfit.updateMask(wsignifMask) + + # Annual slope value for each pixel above threshold. If interested in value for timeperiod, can use .multiply(ee.Image(6)). Can also mask based on offset to limit based on starting vegetation level. + wlfLimit = (wlinearfit.select('scale') # .multiply(ee.Image(6)) + .updateMask(wlinearfit.select('scale').gte(minSlopeWater).Or(wlinearfit.select('scale').lte(ee.Number(0).subtract(minSlopeVeg)))) + ) + + wateranyyearMask = ee.Image(0).where(wanyyear.neq(0), 1) + + wlfLimitanyyearwater = wlfLimit.updateMask(wateranyyearMask) + wlfLimitanyyearwaterLoss = wlfLimitanyyearwater.updateMask(wlfLimitanyyearwater.lt(0)) + wlfLimitanyyearwaterGain = wlfLimitanyyearwater.updateMask(wlfLimitanyyearwater.gt(0)) + + gorwstartmask = gstartmask.select('NDVI').blend(wstartmask.select('NDWI')) + gorwendmask = gendmask.select('NDVI').blend(wendmask.select('NDWI')) + greenorwaterLimitLoss = glfLimitanyyeargreenLoss.blend(wlfLimitanyyearwaterLoss) + greenorwaterLimitGain = glfLimitanyyeargreenGain.blend(wlfLimitanyyearwaterGain) + combinedStartLossGain = gorwstartmask.rename('startgreenwaterIndex').addBands(gorwendmask.rename('endgreenwaterIndex')).addBands( + greenorwaterLimitLoss.rename('lossgreenwaterSlope')).addBands(greenorwaterLimitGain.rename('gaingreenwaterSlope')) + + return combinedStartLossGain + + s2cloudmasked = Albedo().get_masked_s2_collection(ee.Geometry.BBox(*bbox), self.start_date, self.end_date) + vegwatermap = get_map_vegwaterchange(s2cloudmasked) + + data = get_image_collection(ee.ImageCollection(vegwatermap), bbox, self.spatial_resolution, "vegetation water map") + data = data.to_array() + + return data From cc3ec1ef7cca20ba7b75683d1d36ff834956f55d Mon Sep 17 00:00:00 2001 From: weiqi-tori Date: Fri, 25 Oct 2024 14:00:11 +0800 Subject: [PATCH 2/6] update indicator and layer --- city_metrix/layers/__init__.py | 1 + city_metrix/layers/vegetation_water_map.py | 6 ++---- city_metrix/metrics/vegetation_water_change.py | 14 ++++++++++++++ tests/test_layers.py | 5 +++++ 4 files changed, 22 insertions(+), 4 deletions(-) create mode 100644 city_metrix/metrics/vegetation_water_change.py diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py index e701e52b..62538368 100644 --- a/city_metrix/layers/__init__.py +++ b/city_metrix/layers/__init__.py @@ -21,3 +21,4 @@ from .nasa_dem import NasaDEM from .era_5_hottest_day import Era5HottestDay from .impervious_surface import ImperviousSurface +from .vegetation_water_map import VegetationWaterMap diff --git a/city_metrix/layers/vegetation_water_map.py b/city_metrix/layers/vegetation_water_map.py index b936de4e..722e3830 100644 --- a/city_metrix/layers/vegetation_water_map.py +++ b/city_metrix/layers/vegetation_water_map.py @@ -13,7 +13,7 @@ class VegetationWaterMap(Layer): spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) """ - def __init__(self, start_date="2018-01-01", end_date="2022-12-31", spatial_resolution=30, **kwargs): + def __init__(self, start_date="2018-01-01", end_date="2022-12-31", spatial_resolution=10, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date @@ -79,9 +79,6 @@ def sedenomwrap(Im): def eeCdf(t): return ee.Image(0.5).multiply(ee.Image(1).add(ee.Image(t).divide(ee.Image(2).sqrt()).erf())) - def invCdf(p): - return ee.Image(2).sqrt().multiply(ee.Image(p).multiply(2).subtract(1).erfInv()) - # /////green normalized Difference/////// def significanceGreen(gtrend, glinearfit, SampleNumber): mosaicmean = gtrend.mean() @@ -218,6 +215,7 @@ def get_map_vegwaterchange(IC): vegwatermap = get_map_vegwaterchange(s2cloudmasked) data = get_image_collection(ee.ImageCollection(vegwatermap), bbox, self.spatial_resolution, "vegetation water map") + # xarray.dataset to xarray.dataarray data = data.to_array() return data diff --git a/city_metrix/metrics/vegetation_water_change.py b/city_metrix/metrics/vegetation_water_change.py new file mode 100644 index 00000000..54f3ab1c --- /dev/null +++ b/city_metrix/metrics/vegetation_water_change.py @@ -0,0 +1,14 @@ +from geopandas import GeoDataFrame, GeoSeries + +from city_metrix.layers import VegetationWaterMap + + +def vegetation_water_change(zones: GeoDataFrame) -> GeoSeries: + for i in range(len(zones)): + veg_water_map = VegetationWaterMap().get_data(zones.iloc[[i]].total_bounds) + + counts = vegwaterImg.select('startgreenwaterIndex').reduceRegions(collection=boundary,reducer=ee.Reducer.count().setOutputs(['greenorwater2018']),scale=30)#,tileScale=10) + counts = vegwaterImg.select('lossgreenwaterSlope').reduceRegions(collection=counts,reducer=ee.Reducer.count().setOutputs(['greenorwaterLoss']),scale=30)#,tileScale=10) + counts = vegwaterImg.select('gaingreenwaterSlope').reduceRegions(collection=counts,reducer=ee.Reducer.count().setOutputs(['greenorwaterGain']),scale=30)#,tileScale=10) + + return diff --git a/tests/test_layers.py b/tests/test_layers.py index f8e1fe23..41d26ea7 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -24,6 +24,7 @@ TreeCanopyHeight, TreeCover, UrbanLandUse, + VegetationWaterMap, WorldPop ) from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1 @@ -122,6 +123,10 @@ def test_urban_land_use(): data = UrbanLandUse().get_data(BBOX) assert np.size(data) > 0 +def test_vegetation_water_map(): + data = VegetationWaterMap().get_data(BBOX) + assert np.size(data) > 0 + def test_world_pop(): data = WorldPop().get_data(BBOX) assert np.size(data) > 0 From 3e5e916f3ea7fd0402fc908657bfcf020946279b Mon Sep 17 00:00:00 2001 From: weiqi-tori Date: Mon, 11 Nov 2024 14:07:45 +0800 Subject: [PATCH 3/6] update to return one data layer at a time --- city_metrix/layers/vegetation_water_map.py | 26 ++++++++++++------- city_metrix/metrics/__init__.py | 1 + .../metrics/vegetation_water_change.py | 11 ++++---- tests/test_metrics.py | 7 +++++ 4 files changed, 29 insertions(+), 16 deletions(-) diff --git a/city_metrix/layers/vegetation_water_map.py b/city_metrix/layers/vegetation_water_map.py index 722e3830..b89d95eb 100644 --- a/city_metrix/layers/vegetation_water_map.py +++ b/city_metrix/layers/vegetation_water_map.py @@ -10,13 +10,15 @@ class VegetationWaterMap(Layer): Attributes: start_date: starting date for data retrieval end_date: ending date for data retrieval + greenwater_layer: select returned layer from 'startgreenwaterIndex'/'endgreenwaterIndex'/'lossgreenwaterSlope'/'gaingreenwaterSlope' spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) """ - def __init__(self, start_date="2018-01-01", end_date="2022-12-31", spatial_resolution=10, **kwargs): + def __init__(self, start_date="2016-01-01", end_date="2022-12-31", greenwater_layer='startgreenwaterIndex', spatial_resolution=10, **kwargs): super().__init__(**kwargs) self.start_date = start_date self.end_date = end_date + self.greenwater_layer = greenwater_layer self.spatial_resolution = spatial_resolution def get_data(self, bbox): @@ -115,9 +117,9 @@ def senumwrap(Im): wP = ee.Image(1).subtract(eeCdf(wT.abs())) return wP - # function to generate vegetation and water trend and change maps - def get_map_vegwaterchange(IC): + # function to generate vegetation and water trend and change maps + def get_map_vegwaterchange(IC, greenwater_layer): gwic2019 = AnnualIC(IC, '2019') gwic2020 = AnnualIC(IC, '2020') gwic2021 = AnnualIC(IC, '2021') @@ -206,16 +208,20 @@ def get_map_vegwaterchange(IC): gorwendmask = gendmask.select('NDVI').blend(wendmask.select('NDWI')) greenorwaterLimitLoss = glfLimitanyyeargreenLoss.blend(wlfLimitanyyearwaterLoss) greenorwaterLimitGain = glfLimitanyyeargreenGain.blend(wlfLimitanyyearwaterGain) - combinedStartLossGain = gorwstartmask.rename('startgreenwaterIndex').addBands(gorwendmask.rename('endgreenwaterIndex')).addBands( - greenorwaterLimitLoss.rename('lossgreenwaterSlope')).addBands(greenorwaterLimitGain.rename('gaingreenwaterSlope')) - return combinedStartLossGain + if greenwater_layer == 'startgreenwaterIndex': + return gorwstartmask.rename('greenwater_layer') + elif greenwater_layer == 'endgreenwaterIndex': + return gorwendmask.rename('greenwater_layer') + elif greenwater_layer == 'lossgreenwaterSlope': + return greenorwaterLimitLoss.rename('greenwater_layer') + elif greenwater_layer == 'gaingreenwaterSlope': + return greenorwaterLimitGain.rename('greenwater_layer') + s2cloudmasked = Albedo().get_masked_s2_collection(ee.Geometry.BBox(*bbox), self.start_date, self.end_date) - vegwatermap = get_map_vegwaterchange(s2cloudmasked) + vegwatermap = get_map_vegwaterchange(s2cloudmasked, self.greenwater_layer) data = get_image_collection(ee.ImageCollection(vegwatermap), bbox, self.spatial_resolution, "vegetation water map") - # xarray.dataset to xarray.dataarray - data = data.to_array() - return data + return data.greenwater_layer diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index d95cfa5c..b751033f 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -5,3 +5,4 @@ from .urban_open_space import urban_open_space from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing +from .vegetation_water_change import vegetation_water_change diff --git a/city_metrix/metrics/vegetation_water_change.py b/city_metrix/metrics/vegetation_water_change.py index 54f3ab1c..110a3771 100644 --- a/city_metrix/metrics/vegetation_water_change.py +++ b/city_metrix/metrics/vegetation_water_change.py @@ -4,11 +4,10 @@ def vegetation_water_change(zones: GeoDataFrame) -> GeoSeries: - for i in range(len(zones)): - veg_water_map = VegetationWaterMap().get_data(zones.iloc[[i]].total_bounds) + start_counts = VegetationWaterMap(greenwater_layer='startgreenwaterIndex').groupby(zones).count() + loss_counts = VegetationWaterMap(greenwater_layer='lossgreenwaterSlope').groupby(zones).count() + gain_counts = VegetationWaterMap(greenwater_layer='gaingreenwaterSlope').groupby(zones).count() - counts = vegwaterImg.select('startgreenwaterIndex').reduceRegions(collection=boundary,reducer=ee.Reducer.count().setOutputs(['greenorwater2018']),scale=30)#,tileScale=10) - counts = vegwaterImg.select('lossgreenwaterSlope').reduceRegions(collection=counts,reducer=ee.Reducer.count().setOutputs(['greenorwaterLoss']),scale=30)#,tileScale=10) - counts = vegwaterImg.select('gaingreenwaterSlope').reduceRegions(collection=counts,reducer=ee.Reducer.count().setOutputs(['greenorwaterGain']),scale=30)#,tileScale=10) + # TODO: layer generation and zonal stats use different spatial resolutions - return + return (gain_counts - loss_counts) / start_counts diff --git a/tests/test_metrics.py b/tests/test_metrics.py index 8fd42cc3..ae440aaf 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -49,3 +49,10 @@ def test_urban_open_space(): expected_zone_size = ZONES.geometry.size actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size + + +def test_vegetation_water_change(): + indicator = vegetation_water_change(ZONES) + expected_zone_size = ZONES.geometry.size + actual_indicator_size = indicator.size + assert expected_zone_size == actual_indicator_size From d7f1d26f1b3ba73144ff8a4996377f09e8d950d8 Mon Sep 17 00:00:00 2001 From: weiqi-tori Date: Mon, 11 Nov 2024 14:31:29 +0800 Subject: [PATCH 4/6] fix albedo --- city_metrix/layers/albedo.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py index 38983397..f50dd30f 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -13,9 +13,6 @@ class Albedo(Layer): spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) threshold: threshold value for filtering the retrieval """ - S2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED") - S2C = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY") - MAX_CLOUD_PROB = 30 S2_ALBEDO_EQN = '((B*Bw)+(G*Gw)+(R*Rw)+(NIR*NIRw)+(SWIR1*SWIR1w)+(SWIR2*SWIR2w))' @@ -37,8 +34,8 @@ def get_masked_s2_collection(self, roi, start, end): ee.Filter.bounds(roi) )) # .select('B2','B3','B4','B8','B11','B12') - s2 = self.S2.filter(criteria) - s2c = self.S2C.filter(criteria) + s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED").filter(criteria) + s2c = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY").filter(criteria) s2_with_clouds = (ee.Join.saveFirst('cloud_mask').apply(**{ 'primary': ee.ImageCollection(s2), 'secondary': ee.ImageCollection(s2c), From 5f8a297c610130abdf6cd2304db162a9f4e51841 Mon Sep 17 00:00:00 2001 From: weiqi-tori Date: Wed, 18 Dec 2024 16:22:15 +0800 Subject: [PATCH 5/6] rebase on main and return multiple outputs --- city_metrix/layers/albedo.py | 4 ++-- city_metrix/layers/vegetation_water_map.py | 13 ++++++++++--- .../metrics/vegetation_water_change.py | 19 ++++++++++++++++--- tests/test_metrics.py | 18 ++++++++++++++++-- 4 files changed, 44 insertions(+), 10 deletions(-) diff --git a/city_metrix/layers/albedo.py b/city_metrix/layers/albedo.py index f6c2333a..ac6b6b44 100644 --- a/city_metrix/layers/albedo.py +++ b/city_metrix/layers/albedo.py @@ -75,7 +75,7 @@ def _mcc(im): # 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) + .map(self.mask_clouds_and_rescale) ) # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE') s2_with_clouds_ic = ee.ImageCollection(s2_with_clouds) @@ -102,7 +102,7 @@ def calc_s2_albedo(image): 'SWIR2': image.select('B12') } - albedo = image.expression(S2_ALBEDO_EQN, config).double().rename('albedo') + albedo = image.expression(self.S2_ALBEDO_EQN, config).double().rename('albedo') return albedo diff --git a/city_metrix/layers/vegetation_water_map.py b/city_metrix/layers/vegetation_water_map.py index b89d95eb..ba1a5ea5 100644 --- a/city_metrix/layers/vegetation_water_map.py +++ b/city_metrix/layers/vegetation_water_map.py @@ -1,4 +1,3 @@ -import xarray as xr import ee from .layer import Layer, get_image_collection @@ -46,13 +45,16 @@ def AnnualImgGreen(ic, year): bluestThres = bluest.updateMask(bluest.select('NDWI').gte(NDWIthreshold)) greenestnowater = greenest.updateMask(bluestThres.select('NDWI').unmask().Not()) return greenestnowater + def AnnualImgWater(ic, year): bluest = ic.select('NDWI').reduce(ee.Reducer.median()).rename('NDWI').addBands(ee.Image(year).rename('time_start')).float() return bluest + def AnnualImgWatermask(ic, year): bluest = ic.qualityMosaic('NDWI').select('NDWI').addBands(ee.Image(year).rename('time_start')).float() bluestThres = bluest.updateMask(bluest.select('NDWI').gte(NDWIthreshold)) return bluestThres + def AnnualImgGreenmask(ic, year): greenest = ic.qualityMosaic('NDVI').select('NDVI').addBands(ee.Image(year).rename('time_start')).float() greenestThres = greenest.updateMask(greenest.select('NDVI').gte(NDVIthreshold)) @@ -222,6 +224,11 @@ def get_map_vegwaterchange(IC, greenwater_layer): s2cloudmasked = Albedo().get_masked_s2_collection(ee.Geometry.BBox(*bbox), self.start_date, self.end_date) vegwatermap = get_map_vegwaterchange(s2cloudmasked, self.greenwater_layer) - data = get_image_collection(ee.ImageCollection(vegwatermap), bbox, self.spatial_resolution, "vegetation water map") + data = get_image_collection( + ee.ImageCollection(vegwatermap), + bbox, + self.spatial_resolution, + "vegetation water map" + ).greenwater_layer - return data.greenwater_layer + return data diff --git a/city_metrix/metrics/vegetation_water_change.py b/city_metrix/metrics/vegetation_water_change.py index 110a3771..eb4f9e5b 100644 --- a/city_metrix/metrics/vegetation_water_change.py +++ b/city_metrix/metrics/vegetation_water_change.py @@ -2,12 +2,25 @@ from city_metrix.layers import VegetationWaterMap +# TODO: layer generation and zonal stats use different spatial resolutions -def vegetation_water_change(zones: GeoDataFrame) -> GeoSeries: +def vegetation_water_change_gain_area(zones: GeoDataFrame, spatial_resolution=10) -> GeoSeries: + gain_counts = VegetationWaterMap(greenwater_layer='gaingreenwaterSlope').groupby(zones).count() + gain_area = gain_counts * spatial_resolution ** 2 + + return gain_area + + +def vegetation_water_change_loss_area(zones: GeoDataFrame, spatial_resolution=10) -> GeoSeries: + loss_counts = VegetationWaterMap(greenwater_layer='lossgreenwaterSlope').groupby(zones).count() + loss_area = loss_counts * spatial_resolution ** 2 + + return loss_area + + +def vegetation_water_change_gain_loss_ratio(zones: GeoDataFrame) -> GeoSeries: start_counts = VegetationWaterMap(greenwater_layer='startgreenwaterIndex').groupby(zones).count() loss_counts = VegetationWaterMap(greenwater_layer='lossgreenwaterSlope').groupby(zones).count() gain_counts = VegetationWaterMap(greenwater_layer='gaingreenwaterSlope').groupby(zones).count() - # TODO: layer generation and zonal stats use different spatial resolutions - return (gain_counts - loss_counts) / start_counts diff --git a/tests/test_metrics.py b/tests/test_metrics.py index ae440aaf..e3f87830 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -51,8 +51,22 @@ def test_urban_open_space(): assert expected_zone_size == actual_indicator_size -def test_vegetation_water_change(): - indicator = vegetation_water_change(ZONES) +def test_vegetation_water_change_gain_area(): + indicator = vegetation_water_change_gain_area(ZONES) + expected_zone_size = ZONES.geometry.size + actual_indicator_size = indicator.size + assert expected_zone_size == actual_indicator_size + + +def test_vegetation_water_change_loss_area(): + indicator = vegetation_water_change_loss_area(ZONES) + expected_zone_size = ZONES.geometry.size + actual_indicator_size = indicator.size + assert expected_zone_size == actual_indicator_size + + +def test_vegetation_water_change_gain_loss_ratio(): + indicator = vegetation_water_change_gain_loss_ratio(ZONES) expected_zone_size = ZONES.geometry.size actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size From 163abae09822e01f20e7d0e113b0f242821932f3 Mon Sep 17 00:00:00 2001 From: weiqi-tori Date: Wed, 18 Dec 2024 16:30:09 +0800 Subject: [PATCH 6/6] update init file --- city_metrix/metrics/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index b751033f..113e192e 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -5,4 +5,6 @@ from .urban_open_space import urban_open_space from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing -from .vegetation_water_change import vegetation_water_change +from .vegetation_water_change import vegetation_water_change_gain_area +from .vegetation_water_change import vegetation_water_change_loss_area +from .vegetation_water_change import vegetation_water_change_gain_loss_ratio