import argparse
import fsspec
import rasterio as rio # noqa
import xarray as xr
from rasterio.warp import calculate_default_transform
Resampling with Rioxarray (Local storage, NetCDF file, H5NetCDF driver)
def warp_resample(dataset, zoom=0):
from common import earthaccess_args, target_extent
= target_extent[zoom]
te
# Define filepath, driver, and variable information
= earthaccess_args[dataset]
args = f'earthaccess_data/{args["filename"]}'
src # Define source and target projection
= "EPSG:3857"
dstSRS = "EPSG:4326"
srcSRS = height = 256
width
# Specify fsspec caching since default options don't work well for raster data
= {
fsspec_caching "cache_type": "none",
}= fsspec.filesystem("file")
fs with fs.open(src, **fsspec_caching) as f:
# Open dataset
= xr.open_dataset(f, engine="h5netcdf", mask_and_scale=True)[
da "variable"]
args[
]if dataset == "gpm_imerg":
# Transpose and rename dims to align with rioxarray expectations
= da.rename({"lon": "x", "lat": "y"}).transpose("time", "y", "x")
da # Set input dataset projection
= da.rio.write_crs(srcSRS)
da = da.rio.clip_box(
da *te,
=dstSRS,
crs
)# Define affine transformation from input to output dataset
= calculate_default_transform(
dst_transform, w, h
srcSRS,
dstSRS,
da.rio.width,
da.rio.height,*da.rio.bounds(),
=width,
dst_width=height,
dst_height
)# Reproject dataset
return da.rio.reproject(dstSRS, shape=(h, w), transform=dst_transform)
if __name__ == "__main__":
if "get_ipython" in dir():
# Just call warp_resample if running as a Jupyter Notebook
= warp_resample("gpm_imerg")
da else:
# Configure dataset via argpase if running via CLI
= argparse.ArgumentParser(description="Set environment for the script.")
parser
parser.add_argument("--dataset",
="mursst",
defaulthelp="Dataset to resample.",
=["gpm_imerg", "mursst"],
choices
)
parser.add_argument("--zoom",
=0,
defaulthelp="Zoom level for tile extent.",
)= parser.parse_args()
user_args = warp_resample(user_args.dataset, int(user_args.zoom)) da