diff --git a/sample_data/geo_em.d04_LCZ_extent.nc b/sample_data/geo_em.d04_LCZ_extent.nc new file mode 100644 index 0000000..6b4ab7b Binary files /dev/null and b/sample_data/geo_em.d04_LCZ_extent.nc differ diff --git a/sample_data/geo_em.d04_LCZ_params.nc b/sample_data/geo_em.d04_LCZ_params.nc new file mode 100644 index 0000000..cca1211 Binary files /dev/null and b/sample_data/geo_em.d04_LCZ_params.nc differ diff --git a/sample_data/geo_em.d04_NoUrban.nc b/sample_data/geo_em.d04_NoUrban.nc new file mode 100644 index 0000000..e7511ed Binary files /dev/null and b/sample_data/geo_em.d04_NoUrban.nc differ diff --git a/testing/geo_em.d04_LCZ_extent.nc b/testing/geo_em.d04_LCZ_extent.nc index 6b4ab7b..bc3e630 100644 Binary files a/testing/geo_em.d04_LCZ_extent.nc and b/testing/geo_em.d04_LCZ_extent.nc differ diff --git a/testing/geo_em.d04_LCZ_params.nc b/testing/geo_em.d04_LCZ_params.nc index 5f13d72..3b162e7 100644 Binary files a/testing/geo_em.d04_LCZ_params.nc and b/testing/geo_em.d04_LCZ_params.nc differ diff --git a/testing/geo_em.d04_LCZ_params_no_urb_param.nc b/testing/geo_em.d04_LCZ_params_no_urb_param.nc index 0c57cc4..99f405f 100644 Binary files a/testing/geo_em.d04_LCZ_params_no_urb_param.nc and b/testing/geo_em.d04_LCZ_params_no_urb_param.nc differ diff --git a/testing/geo_em.d04_gridinfo.tif b/testing/geo_em.d04_gridinfo.tif deleted file mode 100644 index 252970f..0000000 Binary files a/testing/geo_em.d04_gridinfo.tif and /dev/null differ diff --git a/tests/w2w_test.py b/tests/w2w_test.py index 0c7b38d..6c9354c 100644 --- a/tests/w2w_test.py +++ b/tests/w2w_test.py @@ -23,6 +23,7 @@ from w2w.w2w import _get_lcz_band from w2w.w2w import _get_SW_BW from w2w.w2w import _get_truncated_normal_sample +from w2w.w2w import _get_wrf_grid_info from w2w.w2w import _hgt_resampler from w2w.w2w import _hi_resampler from w2w.w2w import _initialize_urb_param @@ -34,7 +35,6 @@ from w2w.w2w import checks_and_cleaning from w2w.w2w import create_lcz_extent_file from w2w.w2w import create_lcz_params_file -from w2w.w2w import create_wrf_gridinfo from w2w.w2w import expand_land_cat_parents from w2w.w2w import Info from w2w.w2w import main @@ -50,7 +50,6 @@ def _create_info(info_override): 'src_file_clean': '', 'dst_file': '', 'dst_nu_file': '', - 'dst_gridinfo': '', 'dst_lcz_extent_file': '', 'dst_lcz_params_file': '', 'BUILT_LCZ': [], @@ -74,11 +73,10 @@ def test_create_info_dict(): built_lcz=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ) info = Info.from_argparse(args) - # info is Dict, with 9 keys - assert len(info) == 9 + # info is Dict, with 8 keys + assert len(info) == 8 - # Three files are tifs - assert info.dst_gridinfo == 'input/directory/wrf_file_gridinfo.tif' + # Two files are tifs assert info.src_file == 'input/directory/lcz_file.tif' assert info.src_file_clean == 'input/directory/lcz_file_clean.tif' @@ -111,7 +109,7 @@ def test_get_lcz_band_lcz_generator(info_mock, capsys): pytest.param(['--lcz-band', '5'], 5, id='custom band'), ), ) -def test_main_lcz_band_arg(arg, exp, info_mock, capsys, tmpdir): +def test_main_lcz_band_arg(arg, exp, capsys, tmpdir): input_dir = tmpdir.mkdir('input') input_files = ( 'geo_em.d01.nc', @@ -254,7 +252,6 @@ def test_check_lcz_wrf_extent_lcz_too_small(capsys, info_mock): out, _ = capsys.readouterr() assert ('ERROR: LCZ domain should be larger than WRF domain') in out - # TODO maybe add the actual values to check they are correct assert ( 'ERROR: LCZ domain should be larger than WRF domain ' 'in all directions.\nLCZ bounds (xmin, ymin, xmax, ymax): ' @@ -357,24 +354,25 @@ def test_wrf_remove_urban_output_already_exists_is_overwritten(tmpdir, info_mock assert m_time_old != os.path.getmtime(info.dst_nu_file) -def test_create_wrf_gridinfo(tmpdir, info_mock): +def test_get_wrf_grid_info(info_mock): info = info_mock( { - 'dst_nu_file': 'testing/5by5.nc', - 'dst_gridinfo': os.path.join(tmpdir, 'dst_gridinfo.tif'), + 'dst_file': 'sample_data/geo_em.d04.nc', } ) - create_wrf_gridinfo(info) - assert os.listdir(tmpdir) == ['dst_gridinfo.tif'] - tif = rxr.open_rasterio(info.dst_gridinfo) - assert tif.rio.crs.to_proj4() == '+init=epsg:4326' - assert tif.rio.transform() == Affine( - 0.01000213623046875, + wrf_grid_info = _get_wrf_grid_info(info) + crs_string = ( + '+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=5.84999990463257 ' + '+x_0=0 +y_0=0 +R=6370000 +units=m +no_defs' + ) + assert wrf_grid_info['crs'] == crs_string + assert wrf_grid_info['transform'] == Affine( + 1111.7747802734375, 0.0, - -1.4050254821777344, + -838835.115342933, 0.0, - 0.010000228881835938, - 41.480000495910645, + 1111.7747802734375, + 4577176.609447704, ) @@ -437,11 +435,11 @@ def test_get_lcz_arr(info_mock): @pytest.mark.parametrize( ('ucp_key', 'data_sum', 'data_mean'), ( - ('MH_URB2D', 2, 1.1429937), - ('STDH_URB2D', 2, 0.29812142), - ('LB_URB2D', 1, 0.15104416), - ('LF_URB2D', 1, 0.07032267), - ('LP_URB2D', 1, 0.08072148), + ('MH_URB2D', 2, 1.14196181), + ('STDH_URB2D', 2, 0.29785209), + ('LB_URB2D', 1, 0.15090767), + ('LF_URB2D', 1, 0.070259198), + ('LP_URB2D', 1, 0.080648496), ), ) def test_ucp_resampler_output_values_per_paramater( @@ -451,7 +449,8 @@ def test_ucp_resampler_output_values_per_paramater( { 'src_file': 'sample_data/lcz_zaragoza.tif', 'src_file_clean': 'testing/lcz_zaragoza_clean.tif', - 'dst_gridinfo': 'testing/geo_em.d04_gridinfo.tif', + 'dst_file': 'sample_data/geo_em.d04.nc', + 'dst_nu_file': 'testing/geo_em.d04_NoUrban.nc', 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], } ) @@ -476,7 +475,8 @@ def test_ucp_resampler_output_values_per_paramater_frc_threshold(info_mock): { 'src_file': 'sample_data/lcz_zaragoza.tif', 'src_file_clean': 'testing/lcz_zaragoza_clean.tif', - 'dst_gridinfo': 'testing/geo_em.d04_gridinfo.tif', + 'dst_file': 'sample_data/geo_em.d04.nc', + 'dst_nu_file': 'testing/geo_em.d04_NoUrban.nc', 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], } ) @@ -492,7 +492,7 @@ def test_ucp_resampler_output_values_per_paramater_frc_threshold(info_mock): assert ucp_res.shape == (1, 102, 162) # Per parameter, check # max values and non-0 domain average assert np.sum(ucp_res.data == ucp_res.data.max()) == 1 - assert approx(np.mean(ucp_res.data[ucp_res.data > 0])) == 0.428921 + assert approx(np.mean(ucp_res.data[ucp_res.data > 0])) == 0.42892075 # Make sure no nans are present. np.sum(np.isnan(ucp_res.data)) == 0 @@ -503,7 +503,8 @@ def test_hgt_resampler_output_values(info_mock): { 'src_file': 'sample_data/lcz_zaragoza.tif', 'src_file_clean': 'testing/lcz_zaragoza_clean.tif', - 'dst_gridinfo': 'testing/geo_em.d04_gridinfo.tif', + 'dst_file': 'sample_data/geo_em.d04.nc', + 'dst_nu_file': 'testing/geo_em.d04_NoUrban.nc', 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], } ) @@ -520,8 +521,8 @@ def test_hgt_resampler_output_values(info_mock): assert ucp_res.shape == (1, 102, 162) # Check # max values and non-0 domain average - assert np.sum(ucp_res.data == ucp_res.data.max()) == 9 - assert approx(np.mean(ucp_res.data[ucp_res.data > 0])) == 6.7111716 + assert np.sum(ucp_res.data == ucp_res.data.max()) == 10 + assert approx(np.mean(ucp_res.data[ucp_res.data > 0])) == 6.7109852 # Make sure no nans are present. assert np.sum(np.isnan(ucp_res.data)) == 0 @@ -651,7 +652,8 @@ def test_hi_resampler(info_mock): info = info_mock( { 'src_file_clean': 'testing/lcz_zaragoza_clean.tif', - 'dst_gridinfo': 'testing/geo_em.d04_gridinfo.tif', + 'dst_file': 'sample_data/geo_em.d04.nc', + 'dst_nu_file': 'testing/geo_em.d04_NoUrban.nc', 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], } ) @@ -684,15 +686,13 @@ def test_hi_resampler(info_mock): ( ( False, - np.array( - [32.0, 33.0, 35.0, 36.0, 38.0, 39.0, 41.0, 42.0, 43.0, 44.0, 46.0] - ), - np.array([17, 3, 1, 77, 95, 2, 7, 5, 2, 150, 10]), + np.array([32.0, 35.0, 36.0, 38.0, 41.0, 42.0, 43.0, 44.0, 46.0]), + np.array([14, 1, 23, 35, 9, 6, 5, 140, 22]), ), ( True, - np.array([32.0, 33.0, 35.0, 36.0, 38.0, 39.0, 41.0]), - np.array([18, 30, 1, 171, 136, 4, 9]), + np.array([32.0, 33.0, 35.0, 36.0, 38.0, 41.0, np.nan]), + np.array([14, 1, 1, 60, 44, 19, 116]), ), ), ) @@ -706,7 +706,8 @@ def test_lcz_resampler_lcz_nat_mask_on_off_with_lcz15( info = info_mock( { 'src_file_clean': 'testing/lcz_zaragoza_clean.tif', - 'dst_gridinfo': 'testing/geo_em.d04_gridinfo.tif', + 'dst_file': 'sample_data/geo_em.d04.nc', + 'dst_nu_file': 'testing/geo_em.d04_NoUrban.nc', 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15], } ) @@ -722,7 +723,9 @@ def test_lcz_resampler_lcz_nat_mask_on_off_with_lcz15( # With natural masking off, majority filtering also includes # Natural classes lcz_values_def, lcz_counts_def = np.unique(lcz_resampled, return_counts=True) - assert (lcz_values_def == lcz_values).all() + assert ( + lcz_values_def[~np.isnan(lcz_values_def)] == lcz_values[~np.isnan(lcz_values)] + ).all() assert (lcz_counts_def == lcz_counts).all() @@ -733,7 +736,6 @@ def test_add_frc_lu_index_2_wrf(info_mock): 'src_file_clean': 'testing/lcz_zaragoza_clean.tif', 'dst_file': 'sample_data/geo_em.d04.nc', 'dst_nu_file': 'testing/geo_em.d04_NoUrban.nc', - 'dst_gridinfo': 'testing/geo_em.d04_gridinfo.tif', 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], } ) @@ -817,7 +819,6 @@ def test_create_lcz_params_file_attrs_type( input_dir = tmpdir.mkdir('input') shutil.copy(os.path.join('sample_data', 'geo_em.d04.nc'), input_dir) shutil.copy(os.path.join('testing', 'lcz_zaragoza_clean.tif'), input_dir) - shutil.copy(os.path.join('testing', 'geo_em.d04_gridinfo.tif'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_NoUrban.nc'), input_dir) info = info_mock( @@ -825,7 +826,6 @@ def test_create_lcz_params_file_attrs_type( 'src_file_clean': os.path.join(input_dir, 'lcz_zaragoza_clean.tif'), 'dst_file': os.path.join(input_dir, 'geo_em.d04.nc'), 'dst_nu_file': os.path.join(input_dir, 'geo_em.d04_NoUrban.nc'), - 'dst_gridinfo': os.path.join(input_dir, 'geo_em.d04_gridinfo.tif'), 'dst_lcz_params_file': os.path.join(input_dir, 'geo_em.d04_LCZ_params.tif'), 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], } @@ -864,7 +864,6 @@ def test_create_lcz_extent_file(tmpdir, info_mock): 'dst_lcz_extent_file': os.path.join(tmpdir, 'geo_em.d04_LCZ_extent.nc'), 'dst_file': 'sample_data/geo_em.d04.nc', 'dst_nu_file': 'testing/geo_em.d04_NoUrban.nc', - 'dst_gridinfo': 'testing/geo_em.d04_gridinfo.tif', 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], } ) @@ -897,7 +896,7 @@ def test_create_lcz_extent_file(tmpdir, info_mock): assert 'URB_PARAM' not in list(dst_extent.data_vars) # Number of urban pixels within LANDUSEF[12] - assert np.sum(dst_extent['LANDUSEF'].data[0, 12, :, :] == 1) == 369 + assert np.sum(dst_extent['LANDUSEF'].data[0, 12, :, :] == 1) == 255 @pytest.mark.parametrize( @@ -1009,7 +1008,6 @@ def test_checks_and_cleaning_sample_data_all_ok(capsys, tmpdir, info_mock): input_dir = tmpdir.mkdir('input') shutil.copy(os.path.join('sample_data', 'geo_em.d04.nc'), input_dir) shutil.copy(os.path.join('testing', 'lcz_zaragoza_clean.tif'), input_dir) - shutil.copy(os.path.join('testing', 'geo_em.d04_gridinfo.tif'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_NoUrban.nc'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_LCZ_extent.nc'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_LCZ_params.nc'), input_dir) @@ -1019,7 +1017,6 @@ def test_checks_and_cleaning_sample_data_all_ok(capsys, tmpdir, info_mock): 'src_file_clean': os.path.join(input_dir, 'lcz_zaragoza_clean.tif'), 'dst_file': os.path.join(input_dir, 'geo_em.d04.nc'), 'dst_nu_file': os.path.join(input_dir, 'geo_em.d04_NoUrban.nc'), - 'dst_gridinfo': os.path.join(input_dir, 'geo_em.d04_gridinfo.tif'), 'dst_lcz_extent_file': os.path.join(input_dir, 'geo_em.d04_LCZ_extent.nc'), 'dst_lcz_params_file': os.path.join(input_dir, 'geo_em.d04_LCZ_params.nc'), 'BUILT_LCZ': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], @@ -1114,7 +1111,6 @@ def test_checks_and_cleaning_sample_data_check1to5_wrong( input_dir = tmpdir.mkdir('input') shutil.copy(os.path.join('sample_data', 'geo_em.d04.nc'), input_dir) shutil.copy(os.path.join('testing', 'lcz_zaragoza_clean.tif'), input_dir) - shutil.copy(os.path.join('testing', 'geo_em.d04_gridinfo.tif'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_NoUrban_dummy.nc'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_LCZ_extent_dummy.nc'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_LCZ_params_dummy.nc'), input_dir) @@ -1124,7 +1120,6 @@ def test_checks_and_cleaning_sample_data_check1to5_wrong( 'src_file_clean': os.path.join(input_dir, 'lcz_zaragoza_clean.tif'), 'dst_file': os.path.join(input_dir, 'geo_em.d04.nc'), 'dst_nu_file': os.path.join(input_dir, 'geo_em.d04_NoUrban_dummy.nc'), - 'dst_gridinfo': os.path.join(input_dir, 'geo_em.d04_gridinfo.tif'), 'dst_lcz_extent_file': os.path.join( input_dir, 'geo_em.d04_LCZ_extent_dummy.nc' ), @@ -1167,7 +1162,6 @@ def test_checks_and_cleaning_sample_data_check6to9_wrong( input_dir = tmpdir.mkdir('input') shutil.copy(os.path.join('sample_data', 'geo_em.d04.nc'), input_dir) shutil.copy(os.path.join('testing', 'lcz_zaragoza_clean.tif'), input_dir) - shutil.copy(os.path.join('testing', 'geo_em.d04_gridinfo.tif'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_NoUrban_dummy.nc'), input_dir) shutil.copy(os.path.join('testing', 'geo_em.d04_LCZ_extent_dummy.nc'), input_dir) shutil.copy( @@ -1179,7 +1173,6 @@ def test_checks_and_cleaning_sample_data_check6to9_wrong( 'src_file_clean': os.path.join(input_dir, 'lcz_zaragoza_clean.tif'), 'dst_file': os.path.join(input_dir, 'geo_em.d04.nc'), 'dst_nu_file': os.path.join(input_dir, 'geo_em.d04_NoUrban_dummy.nc'), - 'dst_gridinfo': os.path.join(input_dir, 'geo_em.d04_gridinfo.tif'), 'dst_lcz_extent_file': os.path.join( input_dir, 'geo_em.d04_LCZ_extent_dummy.nc' ), diff --git a/w2w/w2w.py b/w2w/w2w.py index 1ec376f..8ad0071 100644 --- a/w2w/w2w.py +++ b/w2w/w2w.py @@ -6,6 +6,7 @@ import rasterio import xarray as xr from rasterio.warp import reproject, Resampling +from rasterio.transform import Affine from scipy.stats import mode, truncnorm import os, sys import argparse @@ -21,7 +22,8 @@ from typing import Union from numpy.typing import NDArray from scipy import stats -from pyproj import CRS +import pyproj +from pyproj import CRS, Transformer if sys.version_info >= (3, 9): # pragma: >=3.9 cover import importlib.metadata as importlib_metadata @@ -152,11 +154,6 @@ def main(argv: Optional[Sequence[str]] = None) -> int: NPIX_NLC=args.NPIX_NLC, ) - print(f'{FBOLD}--> Create temporary WRF grid .tif file ' f'for resampling{FEND}') - create_wrf_gridinfo( - info=info, - ) - print(f'{FBOLD}--> Create LCZ-based geo_em file{FEND}') nbui_max = create_lcz_params_file( info=info, @@ -197,7 +194,6 @@ class Info(NamedTuple): src_file_clean: str dst_file: str dst_nu_file: str - dst_gridinfo: str dst_lcz_extent_file: str dst_lcz_params_file: str BUILT_LCZ: List[int] @@ -215,10 +211,6 @@ def from_argparse(cls, args: argparse.Namespace) -> 'Info': dst_nu_file=os.path.join( args.io_dir, args.wrf_file.replace('.nc', '_NoUrban.nc') ), - # TMP file, will be removed - dst_gridinfo=os.path.join( - args.io_dir, args.wrf_file.replace('.nc', '_gridinfo.tif') - ), dst_lcz_extent_file=os.path.join( args.io_dir, args.wrf_file.replace('.nc', '_LCZ_extent.nc') ), @@ -508,24 +500,87 @@ def wrf_remove_urban( dst_data.to_netcdf(info.dst_nu_file) -# Make WRF grid info available for Resampler (tmp file) -def create_wrf_gridinfo(info: Info) -> None: +# Get WRF grid info for Resampler +# Inspired by: https://fabienmaussion.info/2018/01/06/wrf-projection/ +def _get_wrf_grid_info(info: Info) -> Dict[str, Any]: - # Read gridded WRF data - dst_data = xr.open_dataset(info.dst_nu_file) + # Initialize WGS84 projection + wgs_proj = pyproj.Proj(proj='latlong', datum='WGS84') - # Create simpler WRF grid target. - da_lu = xr.Dataset( - {'LU_INDEX': (['y', 'x'], dst_data['LU_INDEX'][0, :, :].values)}, - coords={ - 'y': dst_data.XLAT_M.values[0, :, 0], - 'x': dst_data.XLONG_M.values[0, 0, :], - }, - ) + # Read gridded WRF data + dst_data = xr.open_dataset(info.dst_file) + + # Projection depends on value of MAP_PROJ + map_proj = int(dst_data.MAP_PROJ) + + # Lambert Conformal Conic + if map_proj == 1: + wrf_proj = pyproj.Proj( + proj='lcc', + units='m', + a=6370000, + b=6370000, + lat_1=dst_data.TRUELAT1, + lat_2=dst_data.TRUELAT2, + lat_0=dst_data.MOAD_CEN_LAT, + lon_0=dst_data.STAND_LON, + ) + # Polar Stereographic + elif map_proj == 2: + hemi = -90.0 if dst_data.TRUELAT1 < 0 else 90.0 + wrf_proj = pyproj.Proj( + proj='stere', + units='m', + a=6370000, + b=6370000, + lat_0=hemi, + lon_0=dst_data.STAND_LON, + lat_ts=dst_data.TRUELAT1, + ) + # Mercator + elif map_proj == 3: + wrf_proj = pyproj.Proj( + proj='merc', + units='m', + a=6370000, + b=6370000, + lon_0=dst_data.STAND_LON, + lat_ts=dst_data.TRUELAT1, + ) + # Latlong - Equidistant Cylindrical + # Follow this: https://github.com/NCAR/wrf-python/blob/ + # 4a9ff241c8f3615b6a5c94e10a945e8a39bdea27/src/wrf/projection.py#L928 + elif map_proj == 6: + wrf_proj = pyproj.Proj( + proj='eqc', + units='m', + a=6370000, + b=6370000, + lon_0=dst_data.STAND_LON, + ) + + # Make transform + transformer_wrf = Transformer.from_proj(wgs_proj, wrf_proj) + e, n = transformer_wrf.transform(dst_data.CEN_LON, dst_data.CEN_LAT) - # Add projection information as attributes, save and read back in. - da_lu.rio.write_crs('epsg:4326', inplace=True) - da_lu.rio.to_raster(info.dst_gridinfo) + # Grid parameters + # https://github.com/fmaussion/salem/blob/ + # d3f2e5e340c2af36c84c82a9de6099c90fba12e8/salem/wrftools.py#L734 + dx, dy = dst_data.DX, dst_data.DY + nx, ny = dst_data.dims['west_east'], dst_data.dims['south_north'] + + # Down left corner of the domain + x0 = -(nx - 1) / 2.0 * dx + e + y0 = -(ny - 1) / 2.0 * dy + n + + wrf_transform = Affine.translation(x0 - dx / 2, y0 - dy / 2) * Affine.scale(dx, dy) + + wrf_grid_info = { + 'crs': wrf_proj.to_proj4(), + 'transform': wrf_transform, + } + + return wrf_grid_info def _get_SW_BW(ucp_table: pd.DataFrame) -> Tuple[pd.Series, pd.Series]: @@ -573,9 +628,12 @@ def _ucp_resampler( '''Helper function to resample lcz ucp data ('FRC_URB2D', 'MH_URB2D', 'STDH_URB2D', 'LB_URB2D', 'LF_URB2D', 'LP_URB2D') to WRF grid''' - # Read gridded data: LCZ and WRF grid + # Read gridded LCZ data src_data = rxr.open_rasterio(info.src_file_clean)[0, :, :] - dst_grid = rxr.open_rasterio(info.dst_gridinfo) + + # Get WRF data and grid info + dst_data = xr.open_dataset(info.dst_nu_file) + wrf_grid_info = _get_wrf_grid_info(info) # Get Street and Building Width SW, BW = _get_SW_BW(ucp_table) @@ -620,11 +678,11 @@ def _ucp_resampler( # Info: https://rasterio.readthedocs.io/en/latest/api/rasterio.warp.html?highlight=reproject(#rasterio.warp.reproject ucp_2_wrf = reproject( lcz_data_da, - dst_grid, + dst_data.LU_INDEX, src_transform=lcz_data_da.rio.transform(), src_crs=lcz_data_da.rio.crs, - dst_transform=dst_grid.rio.transform(), - dst_crs=dst_grid.rio.crs, + dst_transform=wrf_grid_info['transform'], + dst_crs=wrf_grid_info['crs'], resampling=Resampling[RESAMPLE_TYPE], )[0] @@ -647,9 +705,12 @@ def _hgt_resampler( '''Helper function to resample HGT_URB2D (=Area Weighted Mean Building Height ) data to WRF grid''' - # Read gridded data: LCZ and WRF grid + # Read gridded LCZ data src_data = rxr.open_rasterio(info.src_file_clean)[0, :, :] - dst_grid = rxr.open_rasterio(info.dst_gridinfo) + + # Get WRF data and grid info + dst_data = xr.open_dataset(info.dst_nu_file) + wrf_grid_info = _get_wrf_grid_info(info) # Street width extracted from S02012 Building heighht and H2W. SW, BW = _get_SW_BW(ucp_table) @@ -686,22 +747,22 @@ def _hgt_resampler( # Get the aggregated values on WRF grid - nominator ucp_2_wrf_nom = reproject( lcz_data_da_nom, - dst_grid, + dst_data.LU_INDEX, src_transform=lcz_data_da_nom.rio.transform(), src_crs=lcz_data_da_nom.crs, - dst_transform=dst_grid.rio.transform(), - dst_crs=dst_grid.rio.crs, + dst_transform=wrf_grid_info['transform'], + dst_crs=wrf_grid_info['crs'], resampling=Resampling[RESAMPLE_TYPE], )[0].copy() # Get the aggregated values on WRF grid - nominator ucp_2_wrf_denom = reproject( lcz_data_da_denom, - dst_grid, + dst_data.LU_INDEX, src_transform=lcz_data_da_denom.rio.transform(), src_crs=lcz_data_da_denom.crs, - dst_transform=dst_grid.rio.transform(), - dst_crs=dst_grid.rio.crs, + dst_transform=wrf_grid_info['transform'], + dst_crs=wrf_grid_info['crs'], resampling=Resampling[RESAMPLE_TYPE], )[0].copy() @@ -857,9 +918,12 @@ def _hi_resampler( '''Helper function to resample ucp HI_URB2D_URB2D data to WRF grid''' - # Read gridded data: LCZ and WRF grid + # Read gridded LCZ data src_data = rxr.open_rasterio(info.src_file_clean)[0, :, :] - dst_grid = rxr.open_rasterio(info.dst_gridinfo) + + # Get WRF data and grid info + dst_data = xr.open_dataset(info.dst_nu_file) + wrf_grid_info = _get_wrf_grid_info(info) # Get mask of selected built LCZs lcz_arr = _get_lcz_arr(src_data, info) @@ -868,7 +932,7 @@ def _hi_resampler( df_hi = _compute_hi_distribution(info, ucp_table=ucp_table) # Initialize array to store temp values - hi_arr = np.zeros((15, dst_grid.shape[1], dst_grid.shape[2])) + hi_arr = np.zeros((15, dst_data.LU_INDEX.shape[1], dst_data.LU_INDEX.shape[2])) # Loop over the 15 height density classes. for hi_i in range(df_hi.shape[1]): @@ -891,11 +955,11 @@ def _hi_resampler( # Get the aggregated values on WRF grid ucp_2_wrf = reproject( lcz_data_da, - dst_grid, + dst_data.LU_INDEX, src_transform=lcz_data_da.rio.transform(), src_crs=lcz_data_da.rio.crs, - dst_transform=dst_grid.rio.transform(), - dst_crs=dst_grid.rio.crs, + dst_transform=wrf_grid_info['transform'], + dst_crs=wrf_grid_info['crs'], resampling=Resampling[RESAMPLE_TYPE], )[0] @@ -927,10 +991,12 @@ def _lcz_resampler( '''Helper function to resample lcz classes to WRF grid using majority''' - # Read required gridded data, LCZ, WRF grid, and - # original WRF (for original MODIS urban mask) + # Read gridded LCZ data src_data = rxr.open_rasterio(info.src_file_clean)[0, :, :] - dst_grid = rxr.open_rasterio(info.dst_gridinfo) + + # Get WRF data and grid info + dst_data = xr.open_dataset(info.dst_nu_file) + wrf_grid_info = _get_wrf_grid_info(info) # Mask natural LCZs before majority filtering. if LCZ_NAT_MASK: @@ -938,11 +1004,11 @@ def _lcz_resampler( lcz_2_wrf = reproject( src_data, - dst_grid, + dst_data.LU_INDEX, src_transform=src_data.rio.transform(), src_crs=src_data.rio.crs, - dst_transform=dst_grid.rio.transform(), - dst_crs=dst_grid.rio.crs, + dst_transform=wrf_grid_info['transform'], + dst_crs=wrf_grid_info['crs'], resampling=Resampling['mode'], )[0].values @@ -1538,8 +1604,6 @@ def _check_range(darr: NDArray[np.float_], exp_range: List[int]) -> int: ) print('> Cleaning up ... all done!') - if os.path.exists(info.dst_gridinfo): - os.remove(info.dst_gridinfo) if os.path.exists(info.src_file_clean): os.remove(info.src_file_clean)