Skip to content

Commit

Permalink
resolve final conlficts
Browse files Browse the repository at this point in the history
  • Loading branch information
RubenImhoff committed Aug 25, 2021
2 parents efc2d42 + 97352a3 commit 8a9918a
Show file tree
Hide file tree
Showing 2 changed files with 166 additions and 8 deletions.
171 changes: 165 additions & 6 deletions pysteps/io/nwp_importers.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,8 @@
.. autosummary::
:toctree: ../generated/
import_bom_nwp
import_bom_nwp_xr
import_rmi_nwp_xr
import_knmi_nwp
"""

Expand All @@ -89,7 +90,7 @@

@postprocess_import()
def import_bom_nwp_xr(filename, **kwargs):
"""Import a NetCDF with NWP rainfall forecasts regrided to a BoM Rainfields3
"""Import a NetCDF with NWP rainfall forecasts regridded to a BoM Rainfields3
using xarray.
Parameters
Expand All @@ -108,7 +109,7 @@ def import_bom_nwp_xr(filename, **kwargs):

if not NETCDF4_IMPORTED:
raise MissingOptionalDependency(
"netCDF4 package is required to import BoM NWP regrided rainfall "
"netCDF4 package is required to import BoM NWP regridded rainfall "
"products but it is not installed"
)

Expand All @@ -125,7 +126,7 @@ def import_bom_nwp_xr(filename, **kwargs):
# so it needs to be disagregated by time step
varname = kwargs.get("varname", "accum_prcp")
if varname == "accum_prcp":
print("Rainfall values are accumulated. Disagreagating by time step")
print("Rainfall values are accumulated. Disaggregating by time step")
accum_prcp = ds_meta[varname]
precipitation = accum_prcp - accum_prcp.shift({varname_time: 1})
precipitation = precipitation.dropna(varname_time, "all")
Expand Down Expand Up @@ -249,7 +250,165 @@ def _import_bom_nwp_geodata_xr(
"institution": "Commonwealth of Australia, Bureau of Meteorology",
"threshold": _get_threshold_value(da_rainfall.values),
# TODO(_import_bom_rf3_geodata_xr): Remove before final 2.0 version
"yorigin": "lower",
"yorigin": "upper",
"xpixelsize": xpixelsize.values,
"ypixelsize": ypixelsize.values,
}
)

return ds_in


@postprocess_import()
def import_rmi_nwp_xr(filename, **kwargs):
"""Import a NetCDF with NWP rainfall forecasts from RMI using xarray.
Parameters
----------
filename: str
Name of the file to import.
{extra_kwargs_doc}
Returns
-------
out_da : xr.DataArray
A xarray DataArray containing forecast rainfall fields in mm/timestep
imported from a netcdf, and the metadata.
"""

if not NETCDF4_IMPORTED:
raise MissingOptionalDependency(
"netCDF4 package is required to import RMI NWP rainfall "
"products but it is not installed"
)

ds = _import_rmi_nwp_data_xr(filename, **kwargs)
ds_meta = _import_rmi_nwp_geodata_xr(ds, **kwargs)

# rename varname_time (def: time) to t
varname_time = kwargs.get("varname_time", "time")
ds_meta = ds_meta.rename({varname_time: "t"})
varname_time = "t"

# if data variable is named accum_prcp
# it is assumed that NWP rainfall data is accumulated
# so it needs to be disagregated by time step
varname = kwargs.get("varname", "precipitation")
if varname == "accum_prcp":
print("Rainfall values are accumulated. Disaggregating by time step")
accum_prcp = ds_meta[varname]
precipitation = accum_prcp - accum_prcp.shift({varname_time: 1})
precipitation = precipitation.dropna(varname_time, "all")
# update/copy attributes
precipitation.name = "precipitation"
# copy attributes
precipitation.attrs.update({**accum_prcp.attrs})
# update attributes
precipitation.attrs.update({"standard_name": "precipitation_amount"})
else:
precipitation = ds_meta[varname]

return precipitation


def _import_rmi_nwp_data_xr(filename, **kwargs):

varname_time = kwargs.get("varname_time", "time")
chunks = kwargs.get("chunks", {varname_time: 1})

ds_rainfall = xr.open_mfdataset(
filename,
combine="nested",
concat_dim=varname_time,
chunks=chunks,
lock=False,
parallel=True,
)

return ds_rainfall


def _import_rmi_nwp_geodata_xr(
ds_in,
**kwargs,
):

varname = kwargs.get("varname", "precipitation")
varname_time = kwargs.get("varname_time", "time")
projdef = None
if "proj4string" in ds_in.attrs:
projdef = ds_in.proj4string

# get the accumulation period
time = ds_in[varname_time]
# shift the array to calculate the time step
delta_time = time - time.shift({varname_time: 1})
# assuming first valid delta_time is representative of all time steps
time_step = delta_time[1]
time_step = time_step.values.astype("timedelta64[m]")

# get the units of precipitation
units = None
if "units" in ds_in[varname].attrs:
units = ds_in[varname].units
if units in ("kg m-2", "mm"):
units = "mm"
ds_in[varname].attrs.update({"units": units})
# get spatial boundaries and pixelsize
# move to meters if coordinates in kilometers
if "units" in ds_in.x.attrs:
if ds_in.x.units == "km":
ds_in["x"] = ds_in.x * 1000.0
ds_in.x.attrs.update({"units": "m"})
ds_in["y"] = ds_in.y * 1000.0
ds_in.y.attrs.update({"units": "m"})

xmin = ds_in.x.min().values
xmax = ds_in.x.max().values
ymin = ds_in.y.min().values
ymax = ds_in.y.max().values
xpixelsize = abs(ds_in.x[1] - ds_in.x[0])
ypixelsize = abs(ds_in.y[1] - ds_in.y[0])

cartesian_unit = ds_in.x.units

# Add metadata needed by pySTEPS as attrs in X and Y variables

ds_in.x.attrs.update(
{
# TODO: Remove before final 2.0 version
"x1": xmin,
"x2": xmax,
"cartesian_unit": cartesian_unit,
}
)

ds_in.y.attrs.update(
{
# TODO: Remove before final 2.0 version
"y1": ymin,
"y2": ymax,
"cartesian_unit": cartesian_unit,
}
)

# Add metadata needed by pySTEPS as attrs in rainfall variable
da_rainfall = ds_in[varname].isel({varname_time: 0})

ds_in[varname].attrs.update(
{
"transform": None,
"unit": units, # copy 'units' in 'unit' for legacy reasons
"projection": projdef,
"accutime": time_step,
"zr_a": None,
"zr_b": None,
"zerovalue": np.nanmin(da_rainfall),
"institution": "Royal Meteorological Institute of Belgium",
"threshold": _get_threshold_value(da_rainfall.values),
# TODO(_import_bom_rf3_geodata_xr): Remove before final 2.0 version
"yorigin": "upper",
"xpixelsize": xpixelsize.values,
"ypixelsize": ypixelsize.values,
}
Expand Down Expand Up @@ -406,7 +565,7 @@ def _import_knmi_nwp_geodata_xr(

def _get_threshold_value(precip):
"""
Get the the rain/no rain threshold with the same unit, transformation and
Get the rain/no rain threshold with the same unit, transformation and
accutime of the data.
If all the values are NaNs, the returned value is `np.nan`.
Otherwise, np.min(precip[precip > precip.min()]) is returned.
Expand Down
3 changes: 1 addition & 2 deletions pysteps/pystepsrc
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,7 @@
"silent_import": false,
"outputs": {
// path_outputs : path where to save results (figures, forecasts, etc)
"path_outputs": "./",
"path_workdir": "./tmp/"
"path_outputs": "./"
},
"plot": {
// "motion_plot" : "streamplot" or "quiver"
Expand Down

0 comments on commit 8a9918a

Please sign in to comment.