Cif 234 veg water change

city_metrix/layers/__init__.py

@@ -26,3 +26,4 @@
 from .glad_lulc import LandCoverHabitatGlad
 from .glad_lulc import LandCoverHabitatChangeGlad
 from .cams import Cams
+from .vegetation_water_map import VegetationWaterMap

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,8 @@ class Albedo(Layer):
         spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
         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):
         super().__init__(**kwargs)
@@ -20,56 +23,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 = 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),
+            '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 +67,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
-
-

city_metrix/layers/vegetation_water_map.py

@@ -0,0 +1,227 @@
+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
+        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="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):
+        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()))
+
+        # /////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, greenwater_layer):
+            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)
+
+            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, self.greenwater_layer)
+
+        data = get_image_collection(ee.ImageCollection(vegwatermap), bbox, self.spatial_resolution, "vegetation water map")
+
+        return data.greenwater_layer

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

city_metrix/metrics/vegetation_water_change.py

@@ -0,0 +1,13 @@
+from geopandas import GeoDataFrame, GeoSeries
+
+from city_metrix.layers import VegetationWaterMap
+
+
+def vegetation_water_change(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

tests/test_layers.py

@@ -29,6 +29,7 @@
     TreeCanopyHeight,
     TreeCover,
     UrbanLandUse,
+    VegetationWaterMap,
     WorldPop
 )
 from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
@@ -148,6 +149,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