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.

Sign up for GitHub

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

Added treatment of USGS LU classes #99

Merged
merged 8 commits into from
Jan 2, 2023
Merged

Added treatment of USGS LU classes #99

Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
04c93cf
Initial commit
lpilz Oct 31, 2022
d1e8f3a
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Oct 31, 2022
9e5e63d
Fixed tests
lpilz Nov 1, 2022
5b42bfe
Better discernment between USGS and MODIS
lpilz Nov 2, 2022
502c314
Apply suggestions from code review
lpilz Nov 2, 2022
0e75600
Changed confusing naming of is_urban, is_water, is_lake variables
lpilz Nov 16, 2022
5cc43e1
Added test for USGS urbcat removal
lpilz Dec 30, 2022
7010063
Amended LANDUSEF treatment in tests
lpilz Dec 30, 2022
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
84 changes: 52 additions & 32 deletions w2w/w2w.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ def main(argv: Optional[Sequence[str]] = None) -> int:
parser = argparse.ArgumentParser(
description='PURPOSE: Add LCZ-based info to WRF geo_em.dXX.nc\n \n'
'OUTPUT:\n'
'- *_NoUrban.nc: MODIS Urban replaced by surrounding natural LC\n'
'- *_NoUrban.nc: Urban replaced by surrounding natural LC\n'
'- *_LCZ_extent.nc: LCZ urban extent implemented, no LCZ UCPs yet\n'
'- *_LCZ_params.nc: LCZ urban extent + UPC parameter values\n'
'- *_dXX_41.nc: Parent domain files reflecting 41 Land categories',
Expand Down Expand Up @@ -165,9 +165,7 @@ def main(argv: Optional[Sequence[str]] = None) -> int:
LCZ_BAND=LCZ_BAND,
)

print(
f'{FBOLD}--> Replace WRF MODIS urban LC with ' f'surrounding natural LC{FEND}'
)
print(f'{FBOLD}--> Replace WRF urban LC with ' f'surrounding natural LC{FEND}')

wrf_remove_urban(
info=info,
Expand Down Expand Up @@ -431,7 +429,7 @@ def wrf_remove_urban(
NPIX_AREA: int,
) -> None:

'''Remove MODIS urban extent from geo_em*.nc file'''
'''Remove urban extent from geo_em*.nc file'''

# Make a copy of original dst file
dst_data = xr.open_dataset(info.dst_file)
Expand All @@ -443,6 +441,13 @@ def wrf_remove_urban(
lat = dst_data.XLAT_M.squeeze()
lon = dst_data.XLONG_M.squeeze()
orig_num_land_cat = dst_data.NUM_LAND_CAT
is_urban = dst_data.ISURBAN
is_water = dst_data.ISWATER
is_lake = (
dst_data.ISLAKE
if hasattr(dst_data, 'ISLAKE') and dst_data.ISLAKE < orig_num_land_cat
lpilz marked this conversation as resolved.
Show resolved Hide resolved
else None
)

# New arrays to hold data without urban areas
newluse = luse.values.copy()
Expand All @@ -461,14 +466,15 @@ def wrf_remove_urban(
ENDC = '\033[0m'

luf_2D = luf.values.reshape(luf.shape[0], -1)
if orig_num_land_cat > 20: # USING MODIS_LAKE
data_coord['luf_natland'] = list(
(luf_2D[12, :] == 0) & (luf_2D[16, :] == 0) & (luf_2D[20, :] == 0)
data_coord['luf_natland'] = list(
(luf_2D[is_urban - 1, :] == 0) & (luf_2D[is_water - 1, :] == 0)
)
if is_lake is not None:
data_coord['luf_natland'] = data_coord['luf_natland'] & (
luf_2D[is_lake - 1, :] == 0
)
data_coord['luf_urb'] = list((luf_2D[12, :] != 0))
else: # USING MODIS (NO LAKES)
data_coord['luf_natland'] = list((luf_2D[12, :] == 0) & (luf_2D[16, :] == 0))
data_coord['luf_urb'] = list((luf_2D[12, :] != 0))
data_coord['luf_urb'] = list((luf_2D[is_urban - 1, :] != 0))

if NPIX_AREA == None:
NPIX_AREA = NPIX_NLC**2
if NPIX_AREA > luse.size:
Expand All @@ -481,11 +487,14 @@ def wrf_remove_urban(

ikd = using_kdtree(data_coord, min(luse.size, NPIX_AREA))

data_coord['luse_urb'] = data_coord.luse == 13
data_coord['luse_natland'] = (
(data_coord.luse != 13) & (data_coord.luse != 17) & (data_coord.luse != 21)
data_coord['luse_urb'] = data_coord.luse == is_urban
data_coord['luse_natland'] = (data_coord.luse != is_urban) & (
data_coord.luse != is_water
)

if is_lake is not None:
data_coord['luse_natland'] = data_coord['luse_natland'] & (
data_coord.luse != is_lake
)
# Replacing urban pixels with surrounding dominant natural land use category
data_urb = data_coord.where(data_coord.luse_urb).dropna()

Expand Down Expand Up @@ -525,9 +534,9 @@ def wrf_remove_urban(
for iurb_luf in data_luf.index:
i, j = np.unravel_index(iurb_luf, luse.shape)
newluf[int(data_luf.loc[iurb_luf]['newluse']) - 1, i, j] += luf.isel(
south_north=i, west_east=j, land_cat=12
south_north=i, west_east=j, land_cat=is_urban - 1
).values
newluf[12, i, j] = 0.0
newluf[is_urban - 1, i, j] = 0.0

dst_data.LU_INDEX.values[0, :] = newluse[:]
dst_data.LANDUSEF.values[0, :] = newluf[:]
Expand Down Expand Up @@ -1072,12 +1081,13 @@ def _adjust_greenfrac_landusef(

dst_data_orig = xr.open_dataset(info.dst_file)
orig_num_land_cat = dst_data_orig.NUM_LAND_CAT
is_urban = dst_data_orig.ISURBAN

# Adjust GREENFRAC and LANDUSEF
# GREENFRAC is set as average / month from GREENFRAC
# of original MODIS urban pixels
# of original urban pixels
wrf_urb = xr.DataArray(
np.in1d(dst_data_orig['LU_INDEX'][0, :, :].values, [13]).reshape(
np.in1d(dst_data_orig['LU_INDEX'][0, :, :].values, [is_urban]).reshape(
dst_data_orig['LU_INDEX'][0, :, :].shape
),
dims=dst_data_orig['LU_INDEX'][0, :, :].dims,
Expand Down Expand Up @@ -1126,7 +1136,10 @@ def _adjust_greenfrac_landusef(
)
dst_data['LANDUSEF'] = dst_data.LANDUSEF.astype('float32')

luf_attrs['description'] = 'Noah-modified 41-category IGBP-MODIS landuse'
if orig_num_land_cat < 24:
luf_attrs['description'] = 'Noah-modified 41-category ' 'IGBP-MODIS landuse'
else:
luf_attrs['description'] = 'modified 41-category USGS landuse'
for key in luf_attrs.keys():
dst_data['LANDUSEF'].attrs[key] = luf_attrs[key]

Expand Down Expand Up @@ -1157,7 +1170,7 @@ def _add_frc_lu_index_2_wrf(
FRC_THRESHOLD=FRC_THRESHOLD,
)

# Add to geo_em* that that has no MODIS urban
# Add to geo_em* that that has no urban
dst_data = xr.open_dataset(info.dst_nu_file)

# Make a FRC_URB field and store aggregated data.
Expand Down Expand Up @@ -1339,10 +1352,11 @@ def create_lcz_extent_file(info: Info) -> None:

dst_data_orig = xr.open_dataset(info.dst_file)
orig_num_land_cat = dst_data_orig.NUM_LAND_CAT
is_urban = dst_data_orig.ISURBAN
orig_luf_description = dst_data_orig.LANDUSEF.description

lu_index = dst_extent.LU_INDEX.values
lu_index[lu_index >= 31] = 13
lu_index[lu_index >= 31] = is_urban

dst_extent.LU_INDEX.values = lu_index

Expand All @@ -1359,7 +1373,7 @@ def create_lcz_extent_file(info: Info) -> None:
('Time', 'land_cat', 'south_north', 'west_east'),
luf_values[:, :orig_num_land_cat, :, :],
)
dst_extent['LANDUSEF'].values[0, 12, frc_mask] = 1
dst_extent['LANDUSEF'].values[0, is_urban - 1, frc_mask] = 1
dst_extent['LANDUSEF'] = dst_extent.LANDUSEF.astype('float32')

luf_attrs['description'] = orig_luf_description
Expand Down Expand Up @@ -1417,9 +1431,12 @@ def expand_land_cat_parents(info: Info) -> None:
)
da['LANDUSEF'] = da.LANDUSEF.astype('float32')

luf_attrs['description'] = (
'Noah-modified 41-category ' 'IGBP-MODIS landuse'
)
if orig_num_land_cat < 24:
luf_attrs['description'] = (
'Noah-modified 41-category ' 'IGBP-MODIS landuse'
)
else:
luf_attrs['description'] = 'modified 41-category USGS landuse'
for key in luf_attrs.keys():
da['LANDUSEF'].attrs[key] = luf_attrs[key]

Expand Down Expand Up @@ -1456,7 +1473,10 @@ def checks_and_cleaning(info: Info, ucp_table: pd.DataFrame, nbui_max: float) ->
)
ifile = info.dst_nu_file
da = xr.open_dataset(ifile)
if 13 in da.LU_INDEX.values:
dst_data_orig = xr.open_dataset(info.dst_file)
orig_num_land_cat = dst_data_orig.NUM_LAND_CAT
is_urban = dst_data_orig.ISURBAN
if is_urban in da.LU_INDEX.values:
print(
f'{base_text}\n{WARNING} WARNING: Urban land use ' f'still present {ENDC}'
)
Expand All @@ -1469,7 +1489,7 @@ def checks_and_cleaning(info: Info, ucp_table: pd.DataFrame, nbui_max: float) ->
)
ifile = info.dst_lcz_extent_file
da = xr.open_dataset(ifile)
if not 13 in da.LU_INDEX.values:
if not is_urban in da.LU_INDEX.values:
print(
f'{base_text}\n{WARNING} WARNING: LCZ-based urban ' f'extent missing {ENDC}'
)
Expand All @@ -1482,10 +1502,10 @@ def checks_and_cleaning(info: Info, ucp_table: pd.DataFrame, nbui_max: float) ->
)
ifile = info.dst_lcz_params_file
da = xr.open_dataset(ifile)
if 13 in da.LU_INDEX.values:
if is_urban in da.LU_INDEX.values:
print(
f'{base_text}\n{WARNING} WARNING: Urban extent still '
f'defined via LU_INDEX = 13? {ENDC}'
f'defined via LU_INDEX = {is_urban}? {ENDC}'
)
else:
LU_values = np.unique(da.LU_INDEX.values.flatten())
Expand Down Expand Up @@ -1625,7 +1645,7 @@ def _check_range(darr: NDArray[np.float_], exp_range: List[int]) -> int:
)
da_e = xr.open_dataset(info.dst_lcz_extent_file)
da_p = xr.open_dataset(info.dst_lcz_params_file)
da_e_res = xr.where(da_e.LU_INDEX == 13, 1, 0)
da_e_res = xr.where(da_e.LU_INDEX == is_urban, 1, 0)
da_p_res = xr.where(da_p.LU_INDEX >= 31, 1, 0)

if int((da_p_res - da_e_res).sum()) != 0:
Expand Down