Merge pull request #611 from kbarnhart/add_xarray_backends

Add xarray backends for fgmax and fgout gridded formats

examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py

@@ -0,0 +1,87 @@
+import glob
+
+try:
+ import rioxarray
+ import xarray as xr
+except ImportError:
+ "You must install xarray and rioxarray in order to use the xarray backends"
+ raise
+
+from clawpack.geoclaw.xarray_backends import FGMaxBackend, FGOutBackend
+
+# epsg code for lat-lon
+# Optionally, provide an epsg code to assign the associated coordinate system to the file.
+# default behavior assigns no coordinate system. Note that if no coordinate system is provided,
+# it will come in a GIS with coordinates associated with row and column number (not the x and y position
+# encoded in the netcdf).
+
+epsg_code = 4326
+
+# An example of a fgout grid.
+filename = "_output/fgout0001.b0001"
+# provide the .bxxx file if binary format is used or the
+# .qxxx file if ascii format is used.
+# the format, fg number, and frame number are inferred from the filename.
+
+ds = xr.open_dataset(
+ filename,
+ engine=FGOutBackend,
+ backend_kwargs={
+ "epsg": epsg_code,
+ "qmap": "geoclaw",
+ # qmap is the qmap specified to the fgout object in setrun.py see
+ # the following documentation page for more details.
+ # https://www.clawpack.org/dev/fgout.html#specifying-q-out-vars
+ "dry_tolerance": None,
+ # variables that are not eta and B are masked
+ # where h>dry_tolerance. To turn this functionality
+ # off set dry_tolerance = None.
+ },
+)
+# ds is now an xarray object. It can be interacted with directly or written to netcdf using
+ds.to_netcdf("fgout0001_0001.nc")
+
+# It is possible to combine all fgout files into a single netcdf file
+# using xr.open_mfdataset (requires dask) or xr.concat (does not require dask)
+# https://docs.xarray.dev/en/stable/generated/xarray.open_mfdataset.html
+# https://docs.xarray.dev/en/latest/generated/xarray.concat.html
+# for instructions on installing xarray with dask see:
+# https://docs.xarray.dev/en/latest/getting-started-guide/installing.html#instructions
+
+fgout_files = glob.glob("_output/fgout0001.b*")
+
+try:
+ ds_all = xr.open_mfdataset(
+ fgout_files,
+ engine=FGOutBackend,
+ backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"},
+ )
+except ValueError: # if dask is not available, use xr.concat.
+# if dask is not installed xr.open_mfdataset() will fail with something like
+# ValueError: unrecognized chunk manager dask - must be one of: []
+ fgouts = []
+ for filename in fgout_files:
+ ds = xr.open_dataset(
+ filename,
+ engine=FGOutBackend,
+ backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"},
+ )
+ fgouts.append(ds)
+
+ ds_all = xr.concat(fgouts, dim="time")
+
+# save out.
+ds_all.to_netcdf("fgout_all.nc")
+
+# An example of a fgmax grid.
+filename = "_output/fgmax0001.txt"
+ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={"epsg": epsg_code})
+ds.to_netcdf("fgmax0001.nc")
+
+# To see the use of clipping, change the tfinal in setrun to something like 2*3600.0
+# the fgmax0001_clipped.nc will only be the area where the wave arrived within the considered time.
+filename = "_output/fgmax0001.txt"
+ds = xr.open_dataset(
+ filename, engine=FGMaxBackend, backend_kwargs={"epsg": epsg_code, "clip": True}
+)
+ds.to_netcdf("fgmax0001_clipped.nc")