Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Cif 315 improve consistency and standardization of cif code for layers #96

Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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
chrowe marked this conversation as resolved.
Show resolved Hide resolved

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
chrowe marked this conversation as resolved.
Show resolved Hide resolved

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
Loading