Skip to content

Commit

Permalink
Merge pull request #96 from wri/CIF-315-Improve-consistency-and-stand…
Browse files Browse the repository at this point in the history
…ardization-of-CIF-code-for-Layers

Cif 315 improve consistency and standardization of cif code for layers
  • Loading branch information
kcartier-wri authored Dec 17, 2024
2 parents c270fd0 + 0c8b93e commit fbc5a58
Show file tree
Hide file tree
Showing 25 changed files with 476 additions and 248 deletions.
69 changes: 46 additions & 23 deletions city_metrix/layers/albedo.py
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
import ee
import xarray
from dask.diagnostics import ProgressBar

from .layer import Layer, get_utm_zone_epsg, get_image_collection
from .layer import Layer, get_image_collection

class Albedo(Layer):
"""
Expand Down Expand Up @@ -34,19 +32,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(
Expand All @@ -55,23 +67,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):
Expand All @@ -89,20 +107,25 @@ 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_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


23 changes: 14 additions & 9 deletions city_metrix/layers/alos_dsm.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,4 @@
import ee
import xee
import xarray as xr

from .layer import Layer, get_image_collection

Expand All @@ -16,12 +14,19 @@ def __init__(self, spatial_resolution=30, **kwargs):
self.spatial_resolution = spatial_resolution

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
alos_dsm = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2")

alos_dsm_ic = ee.ImageCollection(alos_dsm
.filterBounds(ee.Geometry.BBox(*bbox))
.select('DSM')
.mean()
)

data = get_image_collection(
alos_dsm_ic,
bbox,
self.spatial_resolution,
"ALOS DSM"
).DSM

return data
19 changes: 11 additions & 8 deletions city_metrix/layers/average_net_building_height.py
Original file line number Diff line number Diff line change
@@ -1,9 +1,7 @@
from dask.diagnostics import ProgressBar
import xarray as xr
import xee
import ee

from .layer import Layer, get_utm_zone_epsg, get_image_collection
from .layer import Layer, get_image_collection


class AverageNetBuildingHeight(Layer):
"""
Expand All @@ -23,8 +21,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
18 changes: 11 additions & 7 deletions city_metrix/layers/built_up_height.py
Original file line number Diff line number Diff line change
@@ -1,9 +1,6 @@
from dask.diagnostics import ProgressBar
import xarray as xr
import xee
import ee

from .layer import Layer, get_utm_zone_epsg, get_image_collection
from .layer import Layer, get_image_collection


class BuiltUpHeight(Layer):
Expand All @@ -22,8 +19,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
29 changes: 13 additions & 16 deletions city_metrix/layers/esa_world_cover.py
Original file line number Diff line number Diff line change
@@ -1,9 +1,7 @@
from dask.diagnostics import ProgressBar
from enum import Enum
import xarray as xr
import ee

from .layer import Layer, get_utm_zone_epsg, get_image_collection
from enum import Enum
from .layer import Layer, get_image_collection


class EsaWorldCoverClass(Enum):
Expand Down Expand Up @@ -39,19 +37,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)
Expand Down
14 changes: 9 additions & 5 deletions city_metrix/layers/glad_lulc.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
import xarray as xr
import ee

from .layer import Layer, get_utm_zone_epsg, get_image_collection
from .layer import Layer, get_image_collection


class LandCoverGlad(Layer):
Expand All @@ -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"
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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))
Expand Down
41 changes: 27 additions & 14 deletions city_metrix/layers/high_land_surface_temperature.py
Original file line number Diff line number Diff line change
@@ -1,11 +1,9 @@
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

from shapely.geometry import box
from .land_surface_temperature import LandSurfaceTemperature
from .layer import Layer

class HighLandSurfaceTemperature(Layer):
"""
Expand All @@ -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
Expand All @@ -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):
Expand All @@ -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
Loading

0 comments on commit fbc5a58

Please sign in to comment.