diff --git a/examples/bouss/README.txt b/examples/bouss/README.txt index 047848de2..2e71aa599 100644 --- a/examples/bouss/README.txt +++ b/examples/bouss/README.txt @@ -1,4 +1,6 @@ +GeoClaw Boussinesq solver examples + The radial_flat subdirectory contains one example using the Boussinesq solver introduced in Clawpack v5.10.0. @@ -14,8 +16,12 @@ OpenMP along with MPI. Some flags have to be set as environment variables or directly in the application Makefile, e.g. see the lines commented out in radial_flat/Makefile. +The file setenv.sh illustrates how you might set some environment +variables for the bash shell. + A file petscMPIoptions is also required to set some PETSc parameters for the iterative method used to solve the large sparse linear systems that arise at each refinement level when the Boussinesq equations are solved. One of the environment variables mentioned above points to this file, and a sample is included in this directory. + diff --git a/examples/bouss/radial_flat/README.rst b/examples/bouss/radial_flat/README.rst index 1d76fb818..6c3e6e086 100644 --- a/examples/bouss/radial_flat/README.rst +++ b/examples/bouss/radial_flat/README.rst @@ -8,6 +8,15 @@ A Gaussian hump of water at the origin spreads out over a flat bottom, as specified by the topo file `flat100.txt` (uniform water depth 100 m). The equations are solved in Cartesian coordinates with the SGN equations. +Running the GeoClaw Boussinesq solvers requires PETSc and MPI. +For more details see the documentation + https://www.clawpack.org/bouss2d.html +and + $CLAW/geoclaw/examples/bouss/README.txt +Run + make check +in this directory to check if things seem ok for running this code. + A flagregion specified by `RuledRectangle_Diagonal.data` (that is created by code in `setrun.py` is used to allow flagging for refinement to level 2 only near the diagonal (for abs(x-y) < 1000). The code is set up to diff --git a/examples/tsunami/chile2010/setplot.py b/examples/tsunami/chile2010/setplot.py index 525132748..fc01ec976 100644 --- a/examples/tsunami/chile2010/setplot.py +++ b/examples/tsunami/chile2010/setplot.py @@ -57,18 +57,14 @@ def addgauges(current_data): # Set up for axes in this figure: plotaxes = plotfigure.new_plotaxes('pcolor') - plotaxes.title = 'Surface' - plotaxes.scaled = True - - def fixup(current_data): - import pylab - addgauges(current_data) - t = current_data.t - t = t / 3600. # hours - pylab.title('Surface at %4.2f hours' % t, fontsize=20) - pylab.xticks(fontsize=15) - pylab.yticks(fontsize=15) - plotaxes.afteraxes = fixup + plotaxes.title = 'Surface at time h:m:s' + plotaxes.aspect_latitude = -30. + plotaxes.xticks_fontsize = 10 + plotaxes.yticks_fontsize = 10 + plotaxes.xlabel = 'longitude' + plotaxes.ylabel = 'latitude' + + plotaxes.afteraxes = addgauges # Water plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') @@ -78,6 +74,8 @@ def fixup(current_data): plotitem.pcolor_cmin = -0.2 plotitem.pcolor_cmax = 0.2 plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.8 + plotitem.colorbar_extend = 'both' plotitem.amr_celledges_show = [0,0,0] plotitem.patchedges_show = 1 @@ -114,18 +112,37 @@ def fixup(current_data): # Set up for axes in this figure: plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = 'auto' - plotaxes.ylimits = 'auto' + plotaxes.time_scale = 1/3600. # convert to hours + plotaxes.xlimits = [0, 9] + plotaxes.ylimits = [-0.3, 0.3] plotaxes.title = 'Surface' + plotaxes.title_fontsize = 15 + plotaxes.time_label = 'time (hours)' + plotaxes.ylabel = 'surface elevation (m)' + plotaxes.grid = True + + def add_obs(current_data): + from pylab import plot, legend + gaugeno = current_data.gaugeno + if gaugeno == 32412: + try: + tgauge = TG32412[:,0] / 3600. # convert to hours + plot(tgauge, TG32412[:,1], 'r', linewidth=1.0) + legend(['GeoClaw','Obs'],loc='lower right') + except: pass + + plotaxes.afteraxes = add_obs # Plot surface as blue curve: plotitem = plotaxes.new_plotitem(plot_type='1d_plot') - plotitem.plot_var = 3 + plotitem.plot_var = 3 # eta plotitem.plotstyle = 'b-' + plotitem.kwargs = {'linewidth':1.8} + # Plot topo as green curve: plotitem = plotaxes.new_plotitem(plot_type='1d_plot') - plotitem.show = False + plotitem.show = False # not being used def gaugetopo(current_data): q = current_data.q @@ -137,24 +154,6 @@ def gaugetopo(current_data): plotitem.plot_var = gaugetopo plotitem.plotstyle = 'g-' - def add_zeroline(current_data): - from pylab import plot, legend, xticks, floor, axis, xlabel - t = current_data.t - gaugeno = current_data.gaugeno - - if gaugeno == 32412: - try: - plot(TG32412[:,0], TG32412[:,1], 'r') - legend(['GeoClaw','Obs'],loc='lower right') - except: pass - axis((0,t.max(),-0.3,0.3)) - - plot(t, 0*t, 'k') - n = int(floor(t.max()/3600.) + 2) - xticks([3600*i for i in range(n)], ['%i' % i for i in range(n)]) - xlabel('time (hours)') - - plotaxes.afteraxes = add_zeroline #----------------------------------------- diff --git a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py index 790e8b746..ab47c7588 100644 --- a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py +++ b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py @@ -34,7 +34,6 @@ fgno = 1 # which fgout grid outdir = '_output' -format = 'binary' # format of fgout grid output if 1: # use all fgout frames in outdir: @@ -48,7 +47,8 @@ fgframes = range(1,26) # frames of fgout solution to use in animation # Instantiate object for reading fgout frames: -fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format) +fgout_grid = fgout_tools.FGoutGrid(fgno, outdir) +fgout_grid.read_fgout_grids_data() # Plot one frame of fgout data and define the Artists that will need to diff --git a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py index 4b32c40da..f2d19caad 100644 --- a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py +++ b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py @@ -38,7 +38,6 @@ fgno = 1 # which fgout grid outdir = '_output' -format = 'binary' # format of fgout grid output if 1: # use all fgout frames in outdir: @@ -52,7 +51,8 @@ fgframes = range(1,26) # frames of fgout solution to use in animation # Instantiate object for reading fgout frames: -fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format) +fgout_grid = fgout_tools.FGoutGrid(fgno, outdir) +fgout_grid.read_fgout_grids_data() # Plot one frame of fgout data and define the Artists that will need to diff --git a/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py b/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py index 8b61666d3..a24f5a68d 100644 --- a/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py +++ b/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py @@ -13,6 +13,7 @@ # Instantiate object for reading fgout frames: fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, output_format) +fgout_grid.read_fgout_grids_data() # Plot one frame of fgout data fgframe = 20 diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py new file mode 100644 index 000000000..6d7fa9128 --- /dev/null +++ b/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") diff --git a/src/2d/bouss/amr2.f90 b/src/2d/bouss/amr2.f90 index e8fa9d984..09e2b241b 100644 --- a/src/2d/bouss/amr2.f90 +++ b/src/2d/bouss/amr2.f90 @@ -513,7 +513,7 @@ program amr2 call read_dtopo_settings() ! specifies file with dtopo from earthquake call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_qinit() ! specifies file with dh if this used instead - call set_fgout(rest) ! Fixed grid settings + call set_fgout(rest,nvar) ! Fixed grid settings call setup_variable_friction() ! Variable friction parameter call set_storm() ! Set storm parameters call set_regions() ! Set refinement regions @@ -548,7 +548,7 @@ program amr2 call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_qinit() ! specifies file with dh if this used instead - call set_fgout(rest) ! Fixed grid settings + call set_fgout(rest,nvar) ! Fixed grid settings call setup_variable_friction() ! Variable friction parameter call set_multilayer() ! Set multilayer SWE parameters call set_storm() ! Set storm parameters diff --git a/src/2d/shallow/amr2.f90 b/src/2d/shallow/amr2.f90 index 139e01f1e..7dea01763 100644 --- a/src/2d/shallow/amr2.f90 +++ b/src/2d/shallow/amr2.f90 @@ -491,7 +491,7 @@ program amr2 call read_dtopo_settings() ! specifies file with dtopo from earthquake call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_qinit() ! specifies file with dh if this used instead - call set_fgout(rest) ! Fixed grid settings + call set_fgout(rest,nvar) ! Fixed grid settings call setup_variable_friction() ! Variable friction parameter !call set_multilayer() ! Set multilayer SWE parameters call set_storm() ! Set storm parameters @@ -526,7 +526,7 @@ program amr2 call read_dtopo_settings() ! specifies file with dtopo from earthquake call read_topo_settings(rest) ! specifies topography (bathymetry) files call set_qinit() ! specifies file with dh if this used instead - call set_fgout(rest) ! Fixed grid settings + call set_fgout(rest,nvar) ! Fixed grid settings call setup_variable_friction() ! Variable friction parameter call set_multilayer() ! Set multilayer SWE parameters call set_storm() ! Set storm parameters diff --git a/src/2d/shallow/fgout_module.f90 b/src/2d/shallow/fgout_module.f90 index 981d11ac3..b583df59e 100644 --- a/src/2d/shallow/fgout_module.f90 +++ b/src/2d/shallow/fgout_module.f90 @@ -10,7 +10,7 @@ module fgout_module real(kind=8), pointer :: late(:,:,:) ! Geometry - integer :: num_vars,mx,my,point_style,fgno,output_format,output_style + integer :: mx,my,point_style,fgno,output_format,output_style real(kind=8) :: dx,dy,x_low,x_hi,y_low,y_hi ! Time Tracking and output types @@ -19,11 +19,13 @@ module fgout_module integer, allocatable :: output_frames(:) real(kind=8), allocatable :: output_times(:) - + + integer :: num_vars ! number of variables to interpolate (num_eqn+3) integer :: nqout ! number of q components to output (+1 for eta) - logical, allocatable :: iqout(:) ! which components to output - integer :: bathy_index,eta_index - logical :: dclaw ! False for GeoClaw + logical, allocatable :: q_out_vars(:) ! which components to print + + !logical, allocatable :: iqout(:) ! which components to output + integer :: bathy_index,eta_index,time_index end type fgout_grid @@ -43,14 +45,16 @@ module fgout_module ! Setup routine that reads in the fixed grids data file and sets up the ! appropriate data structures - subroutine set_fgout(rest,fname) + subroutine set_fgout(rest,num_eqn,fname) use amr_module, only: parmunit, tstart_thisrun + use utility_module, only: parse_values implicit none ! Subroutine arguments - logical :: rest ! restart? + logical, intent(in) :: rest ! restart? + integer, intent(in) :: num_eqn character(len=*), optional, intent(in) :: fname ! Local storage @@ -59,6 +63,9 @@ subroutine set_fgout(rest,fname) type(fgout_grid), pointer :: fg real(kind=8) :: ts + integer :: n, iq + real(kind=8) :: values(num_eqn+2) + character(len=80) :: str if (.not.module_setup) then @@ -109,7 +116,25 @@ subroutine set_fgout(rest,fname) read(unit,*) fg%mx, fg%my read(unit,*) fg%x_low, fg%y_low read(unit,*) fg%x_hi, fg%y_hi - read(unit,*) fg%dclaw + + allocate(fg%q_out_vars(num_eqn+2)) ! for q + eta, topo + do iq=1,num_eqn+2 + fg%q_out_vars(iq) = .false. + enddo + read(unit,'(a)') str + call parse_values(str, n, values) + do k=1,n + iq = nint(values(k)) + fg%q_out_vars(iq) = .true. + enddo + write(6,*) '+++ q_out_vars:', fg%q_out_vars + + fg%num_vars = num_eqn + 3 ! also interp topo,eta,time + fg%nqout = 0 ! count how many are to be output + do k=1,num_eqn+2 + if (fg%q_out_vars(k)) fg%nqout = fg%nqout + 1 + enddo + !fg%nqout = fg%nqout + 1 ! for eta ! Initialize next_output_index @@ -193,68 +218,15 @@ subroutine set_fgout(rest,fname) print *,'y grid points my <= 0, set my >= 1' endif - - ! For now, hard-wire with defaults for either GeoClaw or D-Claw - ! need to save q plus topo, eta, t for interp in space-time - - if (fg%dclaw) then - ! For D-Claw: - fg%num_vars = 10 - ! for h, hu, hv, hm, pb, hchi, b_eroded, bathy, eta, time - else - ! GeoClaw: - fg%num_vars = 6 - ! for h, hu, hv, bathy, eta, time - endif - - ! specify which components of q (plus eta?) to output: - ! (eventually this should be set from user input) - - if (fg%num_vars == 6) then - ! GeoClaw - ! indexes used in early and late arrays: - ! 1:3 are q variables, 6 is time - fg%bathy_index = 4 - fg%eta_index = 5 - - allocate(fg%iqout(4)) - fg%iqout(1) = .true. ! output h? - fg%iqout(2) = .true. ! output hu? - fg%iqout(3) = .true. ! output hv? - fg%iqout(4) = .true. ! output eta? - fg%nqout = 0 - do k=1,4 - if (fg%iqout(k)) fg%nqout = fg%nqout + 1 - enddo - endif + fg%bathy_index = num_eqn + 2 + fg%eta_index = num_eqn + 1 + fg%time_index = num_eqn + 3 - if (fg%num_vars == 10) then - ! D-Claw: - ! indexes used in early and late arrays: - ! 1:7 are q variables, 10 is time - fg%bathy_index = 8 - fg%eta_index = 9 - - allocate(fg%iqout(8)) - fg%iqout(1) = .true. ! output h? - fg%iqout(2) = .true. ! output hu? - fg%iqout(3) = .true. ! output hv? - fg%iqout(4) = .true. ! output hm? - fg%iqout(5) = .true. ! output pb? - fg%iqout(6) = .true. ! output hchi? - fg%iqout(7) = .true. ! output beroded? - fg%iqout(8) = .true. ! output eta? - fg%nqout = 0 - do k=1,8 - if (fg%iqout(k)) fg%nqout = fg%nqout + 1 - enddo - endif - write(6,*) '+++ nqout = ',fg%nqout ! Allocate new fixed grid data arrays at early, late time: ! dimension (10,:,:) to work for either GeoClaw or D-Claw - + allocate(fg%early(10, fg%mx,fg%my)) fg%early = nan() @@ -412,7 +384,7 @@ subroutine fgout_interp(fgrid_type,fgrid, & ! save the time this fgout point was computed: - fg_data(num_vars,ifg,jfg) = t + fg_data(fgrid%time_index,ifg,jfg) = t endif ! if fgout point is on this grid @@ -441,7 +413,7 @@ subroutine fgout_write(fgrid,out_time,out_index) ! I/O integer, parameter :: unit = 87 character(len=15) :: fg_filename - character(len=4) :: cfgno, cframeno + character(len=8) :: cfgno, cframeno character(len=8) :: file_format integer :: grid_number,ipos,idigit,out_number,columns integer :: ifg,ifg1, iframe,iframe1 @@ -472,7 +444,7 @@ subroutine fgout_write(fgrid,out_time,out_index) "a10,' format'/,/)" ! Other locals - integer :: i,j,m,iq,k + integer :: i,j,m,iq,k,jq,num_eqn real(kind=8) :: t0,tf,tau, qaug(10) real(kind=8), allocatable :: qeta(:,:,:) real(kind=4), allocatable :: qeta4(:,:,:) @@ -546,76 +518,34 @@ subroutine fgout_write(fgrid,out_time,out_index) endif endif - ! Output the conserved quantities and eta value + ! Output the requested quantities and eta value: + iq = 1 - ! qaug(1:3) are h,hu,hv for both GeoClaw and D-Claw: - if (fgrid%iqout(1)) then - qeta(iq,i,j) = qaug(1) ! h - iq = iq+1 - endif - if (fgrid%iqout(2)) then - qeta(iq,i,j) = qaug(2) ! hu - iq = iq+1 - endif - if (fgrid%iqout(3)) then - qeta(iq,i,j) = qaug(3) ! hv - iq = iq+1 - endif - - if (fgrid%num_vars == 6) then - ! GeoClaw: - if (fgrid%iqout(4)) then - qeta(iq,i,j) = qaug(5) ! eta since qaug(4)=topo - iq = iq+1 - endif - - else if (fgrid%num_vars == 10) then - ! D-Claw: - if (fgrid%iqout(4)) then - qeta(iq,i,j) = qaug(4) ! hm - iq = iq+1 - endif - if (fgrid%iqout(5)) then - qeta(iq,i,j) = qaug(5) ! pb - iq = iq+1 - endif - if (fgrid%iqout(6)) then - qeta(iq,i,j) = qaug(6) ! hchi - iq = iq+1 - endif - if (fgrid%iqout(7)) then - qeta(iq,i,j) = qaug(7) ! b_eroded - iq = iq+1 - endif - if (fgrid%iqout(8)) then - qeta(iq,i,j) = qaug(9) ! eta since qaug(8)=topo + num_eqn = size(fgrid%q_out_vars) - 2 ! last components eta,B + do jq=1,num_eqn+2 + if (fgrid%q_out_vars(jq)) then + qeta(iq,i,j) = qaug(jq) iq = iq+1 endif - endif + enddo + + ! now append eta value: (this is now included in loop above) + !qeta(iq,i,j) = qaug(fgrid%eta_index) + ! NOTE: qaug(fgrid%bathy_index) contains topo B value + enddo enddo - ! Make the file names and open output files - cfgno = '0000' - ifg = fgrid%fgno - ifg1 = ifg - do ipos=4,1,-1 - idigit = mod(ifg1,10) - cfgno(ipos:ipos) = char(ichar('0') + idigit) - ifg1 = ifg1/10 - enddo + ! convert fgrid%fgno to a string of length at least 4 with leading 0's + ! e.g. '0001' or '12345': + write(cfgno, '(I0.4)') fgrid%fgno - cframeno = '0000' - iframe = out_index - iframe1 = iframe - do ipos=4,1,-1 - idigit = mod(iframe1,10) - cframeno(ipos:ipos) = char(ichar('0') + idigit) - iframe1 = iframe1/10 - enddo + ! convert out_index to a string of length at least 4 with leading 0's + write(cframeno, '(I0.4)') out_index - fg_filename = 'fgout' // cfgno // '.q' // cframeno + + fg_filename = 'fgout' // trim(cfgno) // '.q' // trim(cframeno) open(unit,file=fg_filename,status='unknown',form='formatted') @@ -637,14 +567,14 @@ subroutine fgout_write(fgrid,out_time,out_index) if (fgrid%output_format == 3) then ! binary output goes in .b file as full 8-byte (float64): - fg_filename = 'fgout' // cfgno // '.b' // cframeno + fg_filename = 'fgout' // trim(cfgno) // '.b' // trim(cframeno) open(unit=unit, file=fg_filename, status="unknown", & access='stream') write(unit) qeta close(unit) else if (fgrid%output_format == 2) then ! binary output goes in .b file as 4-byte (float32): - fg_filename = 'fgout' // cfgno // '.b' // cframeno + fg_filename = 'fgout' // trim(cfgno) // '.b' // trim(cframeno) open(unit=unit, file=fg_filename, status="unknown", & access='stream') allocate(qeta4(fgrid%nqout, fgrid%mx, fgrid%my)) @@ -671,7 +601,7 @@ subroutine fgout_write(fgrid,out_time,out_index) write(6,*) '*** should be ascii, binary32, or binary64' endif - fg_filename = 'fgout' // cfgno // '.t' // cframeno + fg_filename = 'fgout' // trim(cfgno) // '.t' // trim(cframeno) open(unit=unit, file=fg_filename, status='unknown', form='formatted') ! time, num_eqn+1, num_grids, num_aux, num_dim, num_ghost: write(unit, t_file_format) out_time, fgrid%nqout, 1, 0, 2, 0,file_format @@ -680,10 +610,6 @@ subroutine fgout_write(fgrid,out_time,out_index) print "(a,i4,a,i4,a,e18.8)",'Writing fgout grid #',fgrid%fgno, & ' frame ',out_index,' at time =',out_time - ! Index into qeta for binary output - ! Note that this implicitly assumes that we are outputting only h, hu, hv - ! and will not output more (change num_eqn parameter above) - end subroutine fgout_write diff --git a/src/2d/shallow/surge/storm_module.f90 b/src/2d/shallow/surge/storm_module.f90 index 632cbff90..4c0afe0a4 100644 --- a/src/2d/shallow/surge/storm_module.f90 +++ b/src/2d/shallow/surge/storm_module.f90 @@ -162,7 +162,7 @@ subroutine set_storm(data_file) case(2) rotation => S_rotation case default - stop " *** ERROR *** Roation override invalid." + stop " *** ERROR *** Rotation override invalid." end select read(unit,*) diff --git a/src/2d/shallow/topo_module.f90 b/src/2d/shallow/topo_module.f90 index 8bd1af908..2713e1694 100644 --- a/src/2d/shallow/topo_module.f90 +++ b/src/2d/shallow/topo_module.f90 @@ -439,7 +439,7 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo) integer(kind=8) :: i, j, mtot ! NetCDF Support - character(len=10) :: direction, x_dim_name, x_var_name, y_dim_name, & + character(len=64) :: direction, x_dim_name, x_var_name, y_dim_name, & y_var_name, z_var_name, var_name ! character(len=1) :: axis_string real(kind=8), allocatable :: nc_buffer(:, :), xlocs(:), ylocs(:) @@ -744,8 +744,8 @@ subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy) logical, allocatable :: x_in_dom(:),y_in_dom(:) integer(kind=4) :: dim_ids(2), num_dims, var_type, num_vars, num_dims_tot integer(kind=4), allocatable :: var_ids(:) - character(len=10) :: var_name, x_var_name, y_var_name, z_var_name - character(len=10) :: x_dim_name, y_dim_name + character(len=64) :: var_name, x_var_name, y_var_name, z_var_name + character(len=64) :: x_dim_name, y_dim_name integer(kind=4) :: x_var_id, y_var_id, z_var_id, x_dim_id, y_dim_id verbose = .false. diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py index 83d9fa769..01c132400 100755 --- a/src/python/geoclaw/data.py +++ b/src/python/geoclaw/data.py @@ -9,7 +9,7 @@ - GeoClawData - RefinementData - TopographyData - - FixedGridData + - FGoutData - FGmaxData - DTopoData - QinitData diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py index ed8f675c6..6b957a9d8 100644 --- a/src/python/geoclaw/fgout_tools.py +++ b/src/python/geoclaw/fgout_tools.py @@ -18,7 +18,7 @@ fgout frames, as a single netCDF file - function read_netcdf: Read a netCDF file and return a list of fgout frames, assuming the file contains all the qoi's needed to reconstruct q. -- function read_netcdf_arrays: Read a netCDF file and extract the +- function read_netcdf_arrays: Read a netCDF file and extract the requested quantities of interest as numpy arrays. - print_netcdf_info: Print info about the contents of a netCDF file containing some fgout frames. @@ -52,10 +52,11 @@ def __init__(self, fgout_grid, frameno=None): self._v = None self._s = None self._hss = None + self._hm = None # Define shortcuts to attributes of self.fgout_grid that are the same # for all frames (e.g. X,Y) to avoid storing grid for every frame. - + @property def x(self): return self.fgout_grid.x @@ -63,15 +64,15 @@ def x(self): @property def y(self): return self.fgout_grid.y - + @property def X(self): return self.fgout_grid.X - + @property def Y(self): return self.fgout_grid.Y - + @property def delta(self): return self.fgout_grid.delta @@ -79,15 +80,19 @@ def delta(self): @property def extent_centers(self): return self.fgout_grid.extent_centers - + @property def extent_edges(self): return self.fgout_grid.extent_edges - + @property def drytol(self): return self.fgout_grid.drytol - + + @property + def qmap(self): + return self.fgout_grid.qmap + # Define attributes such as h as @properties with lazy evaluation: # the corresponding array is created and stored only when first # accessed by the user. Those not needed are not created. @@ -96,16 +101,31 @@ def drytol(self): def h(self): """depth""" if self._h is None: - #print('+++ setting _h...') - self._h = self.q[0,:,:] - #print('+++ getting _h...') + q_out_vars = self.fgout_grid.q_out_vars + try: + i_h = q_out_vars.index(self.qmap['h']) + self._h = self.q[i_h,:,:] + except ValueError: # if q_out_vars does not contain qmap['h'], a value error is thrown + try: + i_eta = q_out_vars.index(self.qmap['eta']) + i_B = q_out_vars.index(self.qmap['B']) + self._h = self.q[i_eta,:,:] - self.q[i_B,:,:] + except ValueError: + print('*** Could not find h or eta-B in q_out_vars') + raise AttributeError return self._h @property def hu(self): """momentum h*u""" if self._hu is None: - self._hu = self.q[1,:,:] + q_out_vars = self.fgout_grid.q_out_vars + try: + i_hu = q_out_vars.index(self.qmap['hu']) + self._hu = self.q[i_hu,:,:] + except: + print('*** Could not find hu in q_out_vars') + raise return self._hu @property @@ -121,7 +141,13 @@ def u(self): def hv(self): """momentum h*v""" if self._hv is None: - self._hv = self.q[2,:,:] + q_out_vars = self.fgout_grid.q_out_vars + try: + i_hv = q_out_vars.index(self.qmap['hv']) + self._hv = self.q[i_hv,:,:] + except: + print('*** Could not find hv in q_out_vars') + raise return self._hv @property @@ -137,14 +163,43 @@ def v(self): def eta(self): """surface eta = h+B""" if self._eta is None: - self._eta = self.q[-1,:,:] + q_out_vars = self.fgout_grid.q_out_vars + try: + i_eta = q_out_vars.index(self.qmap['eta']) + self._eta = self.q[i_eta,:,:] + #print('+++qmap["eta"] = %i' % self.qmap["eta"]) + #print('+++i_eta = %i' % i_eta) + except: + try: + i_h = q_out_vars.index(self.qmap['h']) + i_B = q_out_vars.index(self.qmap['B']) + self._eta = self.q[i_h,:,:] + self.q[i_B,:,:] + except: + print('*** Could not find eta or h+B in q_out_vars') + raise return self._eta @property def B(self): """topography""" if self._B is None: - self._B = self.q[-1,:,:] - self.q[0,:,:] + q_out_vars = self.fgout_grid.q_out_vars + try: + i_B = q_out_vars.index(self.qmap['B']) + self._B = self.q[i_B,:,:] + #print('+++qmap["B"] = %i' % self.qmap["B"]) + #print('+++i_B = %i' % i_B) + except: + try: + i_h = q_out_vars.index(self.qmap['h']) + i_eta = q_out_vars.index(self.qmap['eta']) + self._B = self.q[i_eta,:,:] - self.q[i_h,:,:] + #print('+++ computing B: i_h = %i, i_eta = %i' % (i_h,i_eta)) + #print('+++qmap["h"] = %i' % self.qmap["h"]) + #print('+++qmap["eta"] = %i' % self.qmap["eta"]) + except: + print('*** Could not find B or eta-h in q_out_vars') + raise return self._B @property @@ -161,6 +216,86 @@ def hss(self): self._hss = self.h * self.s**2 return self._hss + @property + def huc(self): + """huc - Boussinesq correction to hu""" + if self._huc is None: + q_out_vars = self.fgout_grid.q_out_vars + try: + i_huc = q_out_vars.index(self.qmap['huc']) + self._huc = self.q[i_huc,:,:] + except: + print('*** Could not find huc in q_out_vars') + raise + return self._huc + + @property + def hvc(self): + """hvc - Boussinesq correction to hv""" + if self._hvc is None: + q_out_vars = self.fgout_grid.q_out_vars + try: + i_hvc = q_out_vars.index(self.qmap['hvc']) + self._hvc = self.q[i_hvc,:,:] + except: + print('*** Could not find hvc in q_out_vars') + raise + return self._hvc + + @property + def hm(self): + """dclaw: h * mass fraction""" + if self._hm is None: + q_out_vars = self.fgout_grid.q_out_vars + try: + i_hm = q_out_vars.index(self.qmap['hm']) + self._hm = self.q[i_hm,:,:] + #print('+++qmap["hm"] = %i' % self.qmap["hm"]) + #print('+++i_hm = %i' % i_hm) + except: + print('*** Could not find hm in q_out_vars') + raise + return self._hm + + @property + def pb(self): + """dclaw variable """ + if self._pb is None: + q_out_vars = self.fgout_grid.q_out_vars + try: + i_pb = q_out_vars.index(self.qmap['pb']) + self._pb = self.q[i_pb,:,:] + except: + print('*** Could not find pb in q_out_vars') + raise + return self._pb + + @property + def hchi(self): + """dclaw variable """ + if self._hchi is None: + q_out_vars = self.fgout_grid.q_out_vars + try: + i_hchi = q_out_vars.index(self.qmap['hchi']) + self._hchi = self.q[i_hchi,:,:] + except: + print('*** Could not find hchi in q_out_vars') + raise + return self._hchi + + @property + def bdif(self): + """dclaw variable """ + if self._bdif is None: + q_out_vars = self.fgout_grid.q_out_vars + try: + i_bdif = q_out_vars.index(self.qmap['bdif']) + self._bdif = self.q[i_bdif,:,:] + except: + print('*** Could not find bdif in q_out_vars') + raise + return self._bdif + class FGoutGrid(object): @@ -169,10 +304,28 @@ class FGoutGrid(object): fgout input data and the output generated by a GeoClaw run. """ - def __init__(self,fgno=None,outdir=None,output_format=None): - + def __init__(self,fgno=None,outdir='.',output_format=None, + qmap='geoclaw'): + + + # mapping from variable names to possible values in q_out_vars + if type(qmap) is dict: + self.qmap = qmap + elif qmap == 'geoclaw': + # default for GeoClaw: + self.qmap = {'h':1, 'hu':2, 'hv':3, 'eta':4, 'B':5} + elif qmap == 'geoclaw-bouss': + self.qmap = {'h':1, 'hu':2, 'hv':3, 'huc':4, 'hvc':5, + 'eta':6, 'B':7} + elif qmap == 'dclaw': + self.qmap = {'h':1, 'hu':2, 'hv':3, 'hm':4, 'pb':5, 'hchi':6, + 'bdif':7, 'eta':8, 'B':9} + else: + raise ValueError('Invalid qmap: %s' % qmap) # GeoClaw input values: + # Many of these should be read from fgout_grids.data + # using read_fgout_grids_data before using read_frame self.id = '' # identifier, optional self.point_style = 2 # only option currently supported self.npts = None @@ -184,14 +337,22 @@ def __init__(self,fgno=None,outdir=None,output_format=None): self.nout = None self.fgno = fgno self.outdir = outdir + + # Note output_format will be reset by read_fgout_grids_data: self.output_format = output_format - self.dclaw = False - + + # list of which components to print: + if 'eta' in self.qmap.keys(): + # default: h,hu,hv,eta as in previous versions of GeoClaw + self.q_out_vars = [1,2,3,self.qmap['eta']] + else: + self.q_out_vars = [1,2,3] # if qmap['eta'] not defined + self.drytol = 1e-3 # used for computing u,v from hu,hv # private attributes for those that are only created if # needed by the user: - + self._X = None self._Y = None self._x = None @@ -217,7 +378,7 @@ def y(self): dy = (self.y2 - self.y1)/self.ny self._y = numpy.linspace(self.y1+dy/2, self.y2-dy/2, self.ny) return self._y - + @property def X(self): """2D array X of longitudes (cell centers)""" @@ -239,7 +400,7 @@ def delta(self): dy = (self.y2 - self.y1)/self.ny self._delta = (dx,dy) return self._delta - + @property def extent_centers(self): """Extent of cell centers [xmin,xmax,ymin,ymax]""" @@ -247,7 +408,7 @@ def extent_centers(self): self._extent_centers = [self.x.min(), self.x.max(),\ self.y.min(), self.y.max()] return self._extent_centers - + @property def extent_edges(self): """Extent of cell edges [xmin,xmax,ymin,ymax]""" @@ -256,10 +417,10 @@ def extent_edges(self): self._extent_edges = [self.x.min()-dx/2, self.x.max()+dx/2,\ self.y.min()-dy/2, self.y.max()+dy/2] return self._extent_edges - + # Create plotdata of class clawpack.visclaw.ClawPlotData - # only when needed for reading GeoClaw output, + # only when needed for reading GeoClaw output, # since this must be done after fgno, outdir, output_format # have been specified: @@ -292,7 +453,7 @@ def set_plotdata(self): file_prefix_str = plotdata.file_prefix + '.b' else: file_prefix_str = plotdata.file_prefix + '.q' - + # Contrary to usual default, set save_frames to False so that all # the fgout frames read in are not saved in plotdata.framesoln_dict. # Otherwise this might be a memory hog when making an animation with @@ -312,7 +473,10 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'): self.fgno = fgno assert self.fgno is not None, '*** fgno must be set' - with open(data_file) as filep: + data_path = os.path.join(self.outdir, data_file) + print('Reading fgout grid info from \n %s' % data_path) + + with open(data_path) as filep: lines = filep.readlines() fgout_input = None for lineno,line in enumerate(lines): @@ -324,10 +488,17 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'): if fgout_input is None: raise ValueError('fgout grid fgno = %i not found in %s' \ - % (fgno, data_file)) + % (self.fgno, data_file)) lineno = 0 # next input line - self.output_style = int(fgout_input[lineno].split()[lineno]) + next_value = fgout_input[lineno].split()[lineno] # a string + next_value_int = ('.' not in next_value) # should be True in v5.11 + err_msg = '\n*** Expecting integer output_style next in %s' \ + % data_file \ + + '\n If this is fgout data from before v5.11, try using' \ + + '\n read_fgout_grids_data_pre511' + assert next_value_int, err_msg + self.output_style = int(next_value) lineno += 1 if (self.output_style == 1): # equally spaced times: @@ -354,8 +525,13 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'): lineno += 1 if output_format == 1: self.output_format = 'ascii' + elif output_format == 2: + self.output_format = 'binary32' elif output_format == 3: self.output_format = 'binary' + else: + raise NotImplementedError("fgout not implemented for " \ + + "output_format %i" % output_format) print('Reading input for fgno=%i, point_style = %i ' \ % (self.fgno, self.point_style)) if point_style == 2: @@ -371,7 +547,75 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'): else: raise NotImplementedError("fgout not implemented for point_style %i" \ % point_style) + tokens = fgout_input[lineno].split() + self.q_out_vars = [] + for token in tokens: + try: + self.q_out_vars.append(int(token)) + except: + break + print('Found fgout grid q_out_vars = ',self.q_out_vars) + print('Using this mapping to fgout variable names: ') + print(' qmap = ',self.qmap) + + def read_fgout_grids_data_pre511(self, fgno=None, + data_file='fgout_grids.data'): + """ + For backward compatibility, this reads fgout_grids.data files + in the format used prior to v5.11. + In this case, the following values are used, as set in __init__(): + self.output_style = 1 + self.qmap = 'qmap' + self.q_out_vars = [1,2,3,4] # h,hu,hv,eta + + Read input info for fgout grid number fgno from the data file + fgout_grids.data, which should have been created by setrun.py. + This file now contains info about all fgout grids. + """ + + if fgno is not None: + self.fgno = fgno + assert self.fgno is not None, '*** fgno must be set' + + data_path = os.path.join(self.outdir, data_file) + print('Reading fgout grid info from \n %s' % data_path) + + with open(data_path) as filep: + lines = filep.readlines() + fgout_input = None + for lineno,line in enumerate(lines): + if 'fgno' in line: + if int(line.split()[0]) == self.fgno: + fgout_input = lines[lineno+1:] + #print('Found line %i: %s' % (lineno,line)) + break + + if fgout_input is None: + raise ValueError('fgout grid fgno = %i not found in %s' \ + % (fgno, data_file)) + + self.tstart = float(fgout_input[0].split()[0]) + self.tend = float(fgout_input[1].split()[0]) + self.nout = int(fgout_input[2].split()[0]) + self.point_style = point_style = int(fgout_input[3].split()[0]) + output_format = int(fgout_input[4].split()[0]) + if output_format == 1: + self.output_format = 'ascii' + elif output_format == 3: + self.output_format = 'binary' + print('Reading input for fgno=%i, point_style = %i ' \ + % (self.fgno, self.point_style)) + if point_style == 2: + self.nx = nx = int(fgout_input[5].split()[0]) + self.ny = ny = int(fgout_input[5].split()[1]) + self.x1 = float(fgout_input[6].split()[0]) + self.y1 = float(fgout_input[6].split()[1]) + self.x2 = float(fgout_input[7].split()[0]) + self.y2 = float(fgout_input[7].split()[1]) + else: + raise NotImplementedError("fgout not implemented for point_style %i" \ + % point_style) def write_to_fgout_data(self, fid): """ @@ -387,7 +631,7 @@ def write_to_fgout_data(self, fid): errmsg = 'point_style %s is not implement, only point_style==2' \ % point_style raise NotImplementedError(errmsg) - + if self.output_format == 'ascii': output_format = 1 elif self.output_format == 'binary32': @@ -399,7 +643,7 @@ def write_to_fgout_data(self, fid): raise NotImplementedError(errmsg) - + # write header, independent of point_style: #fid = open(self.input_file_name,'w') fid.write("\n") @@ -407,7 +651,7 @@ def write_to_fgout_data(self, fid): fid.write("%i # output_style\n" \ % self.output_style) - + if self.output_style == 1: assert self.tstart is not None, 'Need to set tstart' assert self.tend is not None, 'Need to set tend' @@ -418,7 +662,7 @@ def write_to_fgout_data(self, fid): elif self.output_style == 2: self.nout = len(self.output_times) fid.write("%i %s # nout\n" % (self.nout, 11*" ")) - + # remove [] and , from list of times: output_times_str = repr(list(self.output_times))[1:-1] output_times_str = output_times_str.replace(',','') @@ -454,8 +698,12 @@ def write_to_fgout_data(self, fid): print(" lower left = (%15.10f,%15.10f)" % (x1,y1)) print(" upper right = (%15.10f,%15.10f)" % (x2,y2)) print(" dx = %15.10e, dy = %15.10e" % (dx,dy)) - - fid.write("%s # dclaw" % self.dclaw) + + # q_out_vars is a list of q components to print, e.g. [1,4] + + format = len(self.q_out_vars) * '%s ' + fid.write(format % tuple(self.q_out_vars)+ " # q_out_vars\n") + fid.write('\n') def read_frame(self, frameno): @@ -465,6 +713,22 @@ def read_frame(self, frameno): from datetime import timedelta + try: + fgoutX = self.X + fgoutY = self.Y + except: + msg = '\n*** Before reading frame, you must set FGoutGrid data,' \ + '\n*** Typically by calling read_fgout_grids_data' + raise ValueError(msg) + + # prior to v5.11, self.read_fgout_grids_data() called here + # rather than raising exception... + print(msg) + print('*** Calling read_fgout_grids_data...') + self.read_fgout_grids_data() + fgoutX = self.X + fgoutY = self.Y + try: fr = self.plotdata.getframe(frameno) except: @@ -480,40 +744,11 @@ def read_frame(self, frameno): fgout_frame.q = state.q fgout_frame.t = state.t - + fgout_frame.frameno = frameno - + X,Y = patch.grid.p_centers[:2] - - try: - fgoutX = self.X - fgoutY = self.Y - except: - # presumably x,y, etc. were not set for this fgout_grid - # (e.g. by read_fgout_grids_data) so infer from this frame and set - - # reset most attributes including private _x etc to None: - self.__init__(fgno=self.fgno,outdir=self.outdir, - output_format=self.output_format) - - print(' Setting grid attributes based on Frame %i:' % frameno) - x = X[:,0] - y = Y[0,:] - dx = x[1] - x[0] - dy = y[1] - y[0] - self.x1 = x[0] - dx/2 - self.x2 = x[-1] + dx/2 - self.y1 = y[0] - dy/2 - self.y2 = y[-1] + dy/2 - self.nx = len(x) - self.ny = len(y) - fgoutX = self.X # will be computed from info above - fgoutY = self.Y # will be computed from info above - print(' Grid edges: ',self.extent_edges) - print(' with %i cells in x and %i cells in y' \ - % (self.nx,self.ny)) - - + if not numpy.allclose(X, fgoutX): errmsg = '*** X read from output does not match fgout_grid.X' raise ValueError(errmsg) @@ -521,7 +756,7 @@ def read_frame(self, frameno): if not numpy.allclose(Y, fgoutY): errmsg = '*** Y read from output does not match fgout_grid.Y' raise ValueError(errmsg) - + return fgout_frame @@ -534,29 +769,29 @@ def make_fgout_fcn_xy(fgout, qoi, method='nearest', """ Create a function that can be called at (x,y) and return the qoi interpolated in space from the fgout array. - - qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to + + qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to an attribute of fgout. - - The function returned takes arguments x,y that can be floats or + + The function returned takes arguments x,y that can be floats or (equal length) 1D arrays of values that lie within the spatial - extent of fgout. - - bounds_error and fill_value determine the behavior if (x,y) is not in + extent of fgout. + + bounds_error and fill_value determine the behavior if (x,y) is not in the bounds of the data, as in scipy.interpolate.RegularGridInterpolator. """ - + from scipy.interpolate import RegularGridInterpolator - + try: q = getattr(fgout,qoi) except: print('*** fgout missing attribute qoi = %s?' % qoi) - + err_msg = '*** q must have same shape as fgout.X\n' \ + 'fgout.X.shape = %s, q.shape = %s' % (fgout.X.shape,q.shape) assert fgout.X.shape == q.shape, err_msg - + x1 = fgout.X[:,0] y1 = fgout.Y[0,:] fgout_fcn1 = RegularGridInterpolator((x1,y1), q, method=method, @@ -584,39 +819,39 @@ def make_fgout_fcn_xyt(fgout1, fgout2, qoi, method_xy='nearest', """ Create a function that can be called at (x,y,t) and return the qoi interpolated in space and time between the two frames fgout1 and fgout2. - - qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to + + qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to an attribute of fgout. - + method_xy is the method used in creating the spatial interpolator, and is passed to make_fgout_fcn_xy. - + method_t is the method used for interpolation in time, currently only 'linear' is supported, which linearly interpolates. - - bounds_error and fill_value determine the behavior if (x,y,t) is not in + + bounds_error and fill_value determine the behavior if (x,y,t) is not in the bounds of the data. - + The function returned takes arguments x,y (floats or equal-length 1D arrays) of values that lie within the spatial extent of fgout1, fgout2 (which are assumed to cover the same uniform grid at different times) - and t should be a float that lies between fgout1.t and fgout2.t. + and t should be a float that lies between fgout1.t and fgout2.t. """ - + assert numpy.allclose(fgout1.X, fgout2.X), \ '*** fgout1 and fgout2 must have same X' assert numpy.allclose(fgout1.Y, fgout2.Y), \ '*** fgout1 and fgout2 must have same Y' - + t1 = fgout1.t t2 = fgout2.t #assert t1 < t2, '*** expected fgout1.t < fgout2.t' - + fgout1_fcn_xy = make_fgout_fcn_xy(fgout1, qoi, method=method_xy, bounds_error=bounds_error, fill_value=fill_value) fgout2_fcn_xy = make_fgout_fcn_xy(fgout2, qoi, method=method_xy, bounds_error=bounds_error, fill_value=fill_value) - + def fgout_fcn(x,y,t): """ Function that can be evaluated at single point or arrays (x,y) @@ -645,15 +880,15 @@ def fgout_fcn(x,y,t): qout = (1-alpha)*qout1 + alpha*qout2 else: raise NotImplementedError('method_t = %s not supported' % method_t) - + return qout - + return fgout_fcn # =============================== # Functions for writing a set of fgout frames as a netCDF file, and -# reading such a file: - +# reading such a file: + def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc', qois = ['h','hu','hv','eta'], datatype='f4', include_B0=False, include_Bfinal=False, @@ -662,63 +897,63 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc', Write a list of fgout frames (at different times on the same rectangular grid) to a single netCDF file, with some metadata and the topography, if desired. - + fgout_frames should be a list of FGoutFrame objects, all of the same size and at increasing times. - + fname_nc is the name of the file to write. - + qois is a list of strings, the quantities of interest to include in file. - This could include any of: + This could include any of: 'h', 'eta', 'hu', 'hv', 'u', 'v', 's', 'hss', 'B'. All other quantities can be computed from h, hu, hv, eta, the original fgout variables from GeoClaw, but for some applications you might only want to save 'h' and 's', for example. - + datatype should be 'f4' [default] or 'f8', specifying bytes per qoi value. 'f8' has full precision of the original data, but the file will be twice as large and may not be needed for downstream applications. - + Note that the topography B = eta - h, so it is not necessary to store all three of these. Also, B is often the same for all frames, so rather than storing B at each frame as a qoi, two other options are also provided (and then storing eta or h for all frames allows calculating the other): - + include_Bfinal: If True, include the topography B array from the final frame as the Bfinal array. - + include_B0: If True, include the topography B array from the first frame - as the B0 array. This is only useful if, e.g., the first frame is initial + as the B0 array. This is only useful if, e.g., the first frame is initial topography before co-seismic deformation, and at later times the topography - is always equal to Bfinal. - + is always equal to Bfinal. + `description` is a string that will be added as metadata. - A metadata field `history` will also be added, which includes the - time the file was created and the path to the directory where it was made. - """ - + A metadata field `history` will also be added, which includes the + time the file was created and the path to the directory where it was made. + """ + import netCDF4 from datetime import datetime, timedelta from cftime import num2date, date2num import time timestr = time.ctime(time.time()) # current time for metadata - + fg_times = numpy.array([fg.t for fg in fgout_frames]) - + if verbose: print('Creating %s with fgout frames at times: ' % fname_nc) print(fg_times) - + fg0 = fgout_frames[0] x = fg0.x y = fg0.y - + xs = numpy.array([fg.x for fg in fgout_frames]) ys = numpy.array([fg.y for fg in fgout_frames]) # assert same for all times - + units = {'h':'meters', 'eta':'meters', 'hu':'m^2/s', 'hv':'m^2/s', 'u':'m/s', 'v':'m/s', 's':'m/s', 'hss':'m^3/s^2', 'B':'meters'} @@ -737,12 +972,12 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc', latitudes = rootgrp.createVariable('lat','f8',('lat',)) latitudes[:] = y latitudes.units = 'degrees_north' - + time = rootgrp.createDimension('time', len(fg_times)) times = rootgrp.createVariable('time','f8',('time',)) times[:] = fg_times times.units = 'seconds' - + if 0: # Could make times be datetimes relative to some event time, e.g.: times.units = 'seconds since 1700-01-26 21:00:00.0' @@ -750,7 +985,7 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc', dates = [datetime(1700,1,26,21) + timedelta(seconds=ss) \ for ss in fg_times] times[:] = date2num(dates,units=times.units,calendar=times.calendar) - + if include_B0: B0 = rootgrp.createVariable('B0',datatype,('lon','lat',)) B0[:,:] = fg0.B @@ -761,7 +996,7 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc', Bfinal = rootgrp.createVariable('Bfinal',datatype,('lon','lat',)) Bfinal[:,:] = fg_final.B Bfinal.units = 'meters' - + for qoi in qois: qoi_frames = [getattr(fgout,qoi) for fgout in fgout_frames] qoi_var = rootgrp.createVariable(qoi,datatype,('time','lon','lat',)) @@ -770,7 +1005,7 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc', def get_as_array(var, rootgrp, verbose=True): """ - Utility function to retrieve variable from netCDF file and convert to + Utility function to retrieve variable from netCDF file and convert to numpy array. """ a = rootgrp.variables.get(var, None) @@ -779,7 +1014,7 @@ def get_as_array(var, rootgrp, verbose=True): return numpy.array(a) else: if verbose: print(' Did not find %s' % var) - return None + return None def read_netcdf_arrays(fname_nc, qois, verbose=True): @@ -790,19 +1025,19 @@ def read_netcdf_arrays(fname_nc, qois, verbose=True): 'h', 'eta', 'hu', 'hv', 'u', 'v', 's', 'hss', 'B'. qois can also include the time-independent 'B0' and/or 'Bfinal' - - Returns + + Returns x, y, t, qoi_arrays where x,y define the longitude, latitudes, t is the times of the frames, and qoi_arrays is a dictionary indexed by the strings from qois. - + Each dict element is an array with shape (len(t), len(x), len(y)) for time-dependent qoi's, or (len(x), len(y)) for B0 or Bfinal, or None if that qoi was not found in the netCDF file. """ - + import netCDF4 with netCDF4.Dataset(fname_nc, 'r') as rootgrp: @@ -814,17 +1049,17 @@ def read_netcdf_arrays(fname_nc, qois, verbose=True): x = get_as_array('lon', rootgrp, verbose) y = get_as_array('lat', rootgrp, verbose) t = get_as_array('time', rootgrp, verbose) - + vars = list(rootgrp.variables) - + qoi_arrays = {} for qoi in qois: qoi_array = get_as_array(qoi, rootgrp, verbose) qoi_arrays[qoi] = qoi_array - + return x,y,t,qoi_arrays - - + + def read_netcdf(fname_nc, fgout_grid=None, verbose=True): """ @@ -833,7 +1068,7 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True): the qoi's 'h','hu','hv','eta' required to reconstruct the q array as output by GeoClaw. """ - + import netCDF4 with netCDF4.Dataset(fname_nc, 'r') as rootgrp: @@ -847,7 +1082,7 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True): t = get_as_array('time', rootgrp, verbose) vars = list(rootgrp.variables) - + for qoi in ['h','hu','hv','eta']: errmsg = '*** Cannot reconstruct fgout frame without %s' % qoi @@ -856,14 +1091,14 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True): hu = get_as_array('hu', rootgrp, verbose) hv = get_as_array('hv', rootgrp, verbose) eta = get_as_array('eta', rootgrp, verbose) - + if (x is None) or (y is None): print('*** Could not create grid') else: X,Y = numpy.meshgrid(x,y) - + fgout_frames = [] - + for k in range(eta.shape[0]): fgout = FGoutFrame(fgout_grid=fgout_grid, frameno=k) fgout.x = x @@ -877,9 +1112,9 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True): fgout.X = X fgout.Y = Y fgout_frames.append(fgout) - + print('Created fgout_frames as list of length %i' % len(fgout_frames)) - + return fgout_frames def print_netcdf_info(fname_nc): @@ -895,7 +1130,7 @@ def print_netcdf_info(fname_nc): t = get_as_array('time', rootgrp, verbose=False) vars = list(rootgrp.variables) - + print('===================================================') print('netCDF file %s contains:' % fname_nc) print('description: \n', rootgrp.description) @@ -905,5 +1140,3 @@ def print_netcdf_info(fname_nc): print('%i times from %.3f to %.3f' % (len(t),t[0],t[-1])) print('variables: ',vars) print('===================================================') - - diff --git a/src/python/geoclaw/surge/plot.py b/src/python/geoclaw/surge/plot.py index 0464eebac..074fc937c 100644 --- a/src/python/geoclaw/surge/plot.py +++ b/src/python/geoclaw/surge/plot.py @@ -440,7 +440,7 @@ def add_pressure(plotaxes, bounds=None, plot_type='pcolor', shrink=1.0): pass -def add_land(plotaxes, plot_type='pcolor', bounds=None): +def add_land(plotaxes, plot_type='pcolor', bounds=[0, 50]): """Add plotitem for land""" if plot_type == 'pcolor': @@ -448,9 +448,8 @@ def add_land(plotaxes, plot_type='pcolor', bounds=None): plotitem.show = True plotitem.plot_var = geoplot.land plotitem.pcolor_cmap = land_cmap - if bounds is not None: - plotitem.pcolor_cmin = bounds[0] - plotitem.pcolor_cmax = bounds[1] + plotitem.pcolor_cmin = bounds[0] + plotitem.pcolor_cmax = bounds[1] plotitem.add_colorbar = False plotitem.amr_celledges_show = [0] * 10 plotitem.amr_patchedges_show = [1, 1, 1, 1, 1, 0, 0] diff --git a/src/python/geoclaw/surge/storm.py b/src/python/geoclaw/surge/storm.py index 2eab6f110..5b0992ee4 100644 --- a/src/python/geoclaw/surge/storm.py +++ b/src/python/geoclaw/surge/storm.py @@ -1169,7 +1169,7 @@ def plot(self, ax, *args, t_range=None, categorization=None, coordinate_system=2, fill_alpha=0.25, fill_color='red', **kwargs): """Plot this storm's track in the given axes object - + :Input: - *ax* (matplotlib.pyplot.axes) Axes to plot into. - *t_range* (list) Time range to plot the track for. If None then use @@ -1286,7 +1286,6 @@ def plot(self, ax, *args, t_range=None, categorization=None, # ========================================================================= # Other Useful Routines - def category(self, categorization="NHC", cat_names=False): r"""Categorizes storm based on relevant storm data @@ -1339,9 +1338,7 @@ def category(self, categorization="NHC", cat_names=False): 12: "Hurricane"} elif categorization.upper() == "NHC": - # TODO: Change these to m/s (knots are how these are defined). - # Definitely not in the correct format now - # TODO: Add TD and TS designations + # NHC uses knots speeds = units.convert(self.max_wind_speed, "m/s", "knots") category = (np.zeros(speeds.shape) + (speeds < 30) * -1 + diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py new file mode 100644 index 000000000..326d0d908 --- /dev/null +++ b/src/python/geoclaw/xarray_backends.py @@ -0,0 +1,583 @@ +r""" +xarray backends module: $CLAW/geoclaw/src/python/geoclaw/xarray_backends.py + +Xarray backends for GeoClaw fixed grids and fgmax grids. + +These only work for point_style = 2 (uniform regular) and have a dependency on +xarray and rioxarray. Xarray provides the core datastructure and rioxarray +provides an interface to rasterio, used to assign geospatial projection +information. + +- https://docs.xarray.dev/en/stable/index.html +- https://corteva.github.io/rioxarray/stable/readme.html +- https://rasterio.readthedocs.io/en/latest/index.html + +The expectation is that you run commands that use these backends from the same +directory you run "make output". + +Includes: + +- class FGMaxBackend: Xarray backend for fgmax grids. +- class FGOutBackend: Xarray backend for fgout grids. + +Usage: + +.. code-block:: python + import rioxarray # activate the rio accessor + import xarray as xr + from clawpack.geoclaw.xarray_backends import FGOutBackend, FGMaxBackend + + # 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={ + 'qmap': 'geoclaw', + 'epsg':epsg_code + }) + # ds is now an xarray object. It can be interacted with directly or written to netcdf using + ds.write_netcdf('filename.nc') + + # A 'qmap' backend_kwargs is required as it indicates the qmap used for + # selecting elements of q to include in fgout. See + # https://www.clawpack.org/dev/fgout.html#specifying-q-out-vars + + # Optionally, provide an epsg code to assign the associated coordinate system to the file. + # default behavior assigns no coordinate system. + + # An example of a fgmax grid. + filename = "_output/fgmax0001.txt" + ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={'epsg':epsg_code}) + + +Dimensions: + +Files opened with FGOutBackend will have dimensions (time, y, x). +Files opened with FGMaxBackend will have dimensions (y, x). + +Variable naming: + +For fixed grid files, the dataset will have the elements of q specified +by qmap and q_out_vars. This may be as many as the full list described below or +as few as specified by the user. + +For fixed grid geoclaw files, the full set of variables is: +- h +- hu +- hv +- eta +- B + +For fixed grid geoclaw-bouss files, the full set of variables is: +- h +- hu +- hv +- huc +- hvc +- eta +- B + +Fixed grid dclaw files, the full set of variables is: +- h +- hu +- hv +- hm +- pb +- hchi +- bdif +- eta +- B + +Depending on the number of variables specified in the setrun.py fgmax files will +have a portion of the following variables: + +If rundata.fgmax_data.num_fgmax_val == 1 + +- arrival_time, Wave arrival time (based on eta>sea_level + fg.arrival_tol) +- h_max, Maximum water depth +- eta_max, Maximum water surface elevation +- h_max_time, Time of maximum water depth +- B, Basal topography at the first time fgmax first monitored maximum amr level +- level, Maximum amr level + +If rundata.fgmax_data.num_fgmax_val == 2: + +- s_max, Maximum velocity +- s_max_time, Time of maximum velocity + +If rundata.fgmax_data.num_fgmax_val == 5: + +- hs_max, Maximum momentum +- hs_max_time, Time of maximum momentum +- hss_max, Maximum momentum flux +- hss_max_time, Time of maximum momentum flux +- h_min, Minimum depth +- h_min_time, Time of minimum depth + + +See the following links for additional information about xarray Backends. + +- https://docs.xarray.dev/en/stable/generated/xarray.backends.BackendEntrypoint.html#xarray.backends.BackendEntrypoint +- https://docs.xarray.dev/en/stable/generated/xarray.open_dataset.html +""" + +import os + +import numpy as np + +try: + import rioxarray # activate the rio accessor + import xarray as xr +except ImportError: + raise ImportError( + "rioxarray and xarray are required to use the FGOutBackend and FGMaxBackend" + ) + +from clawpack.geoclaw import fgmax_tools, fgout_tools +from xarray.backends import BackendEntrypoint + +_qunits = { + "h": "meters", + "hu": "meters squared per second", + "hv": "meters squared per second", + "hm": "meters", + "pb": "newton per meter squared", + "hchi": "meters", + "bdif": "meters", + "eta": "meters", + "B": "meters", + "huc": "varies", + "hvc": "varies", +} + +time_unit = "seconds" +reference_time = "model start" +space_unit = "meters" +nodata = np.nan + + +def _prepare_var(data, mask, dims, units, nodata, long_name): + "Reorient array into the format xarray expects and generate the mapping." + if mask is not None: + data[mask] = nodata + data = data.T + data = np.flipud(data) + mapping = ( + dims, + data, + {"units": units, "_FillValue": nodata, "long_name": long_name}, + ) + return mapping + + +class FGOutBackend(BackendEntrypoint): + "Xarray Backend for Clawpack fixed grid format." + + def open_dataset( + self, + filename, # path to fgout file. + qmap="geoclaw", # qmap value for FGoutGrid ('geoclaw', 'dclaw', or 'geoclaw-bouss') + epsg=None, # epsg code + dry_tolerance=0.001, + # dry tolerance used for for masking elements of q that are not + # eta or B. Default behavior is to mask all elements of q except for + # eta and B where h>0.001. + # used only if h or eta-B is available based on q_out_vars. + # if dry_tolerance = None, no masking is applied to any variable. + drop_variables=None, # name of any elements of q to drop. + ): + + if drop_variables is None: + drop_variables = [] + + full_path = os.path.abspath(filename) + filename = os.path.basename(full_path) + outdir = os.path.basename(os.path.dirname(full_path)) + + # filename has the format fgoutXXXX.qYYYY (ascii) + # or fgoutXXXX.bYYYY (binary) + # where XXXX is the fixed grid number and YYYY is the frame + # number. + type_code = filename.split(".")[-1][0] + fgno = int(filename.split(".")[0][-4:]) + frameno = int(filename.split(".")[-1][-4:]) + if type_code == "q": # TODO, is this correct? + output_format = "ascii" + elif type_code == "b": + output_format = "binary32" # format of fgout grid output + else: + raise ValueError("Invalid FGout output format. Must be ascii or binary.") + + fgout_grid = fgout_tools.FGoutGrid( + fgno=fgno, outdir=outdir, output_format=output_format, qmap=qmap + ) + + if fgout_grid.point_style != 2: + raise ValueError("FGOutBackend only works with fg.point_style=2") + + fgout_grid.read_fgout_grids_data() + fgout = fgout_grid.read_frame(frameno) + + time = fgout.t + # both come in ascending. flip to give expected order. + x = fgout.x + y = np.flipud(fgout.y) + nj = len(x) + ni = len(y) + + # mask based on dry tolerance + if dry_tolerance is not None: + try: + mask = fgout.h < dry_tolerance + # internally fgout_tools will try fgou.eta-fgout.B if h is not present. + except AttributeError: + print("FGOutBackend: No h, eta, or B. No mask applied.") + mask = np.ones((nj, ni), dtype=bool) + + # create data_vars dictionary + data_vars = {} + for i, i_var in enumerate(fgout_grid.q_out_vars): + + # construct variable + + # Find the varname in fgout.qmap associated with + # q_out_vars[i] + varname = None + for name, index in fgout_grid.qmap.items(): + if i_var == index: + varname = name + + # construct xarray dataset if varname not in drop vars. + if varname not in drop_variables: + + Q = fgout.q[i, :, :] + + # mask all but eta based on h presence. + if (varname not in ("B", "eta")) and dry_tolerance is not None: + Q[mask] = nodata + + # to keep xarray happy, need to transpose and flip ud. + Q = Q.T + Q = np.flipud(Q) + Q = Q.reshape((1, ni, nj)) # reshape to add a time dimension. + + data_array_attrs = {"units": _qunits[varname], "_FillValue": nodata} + + data_vars[varname] = ( + [ + "time", + "y", + "x", + ], + Q, + data_array_attrs, + ) + + ds_attrs = {"description": "Clawpack model output"} + + ds = xr.Dataset( + data_vars=data_vars, + coords=dict( + x=(["x"], x, {"units": space_unit}), + y=(["y"], y, {"units": space_unit}), + time=("time", [time], {"units": "seconds"}), + reference_time=reference_time, + ), + attrs=ds_attrs, + ) + + if epsg is not None: + ds.rio.write_crs( + epsg, + inplace=True, + ).rio.set_spatial_dims( + x_dim="x", + y_dim="y", + inplace=True, + ).rio.write_coordinate_system(inplace=True) + # https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Spatial-dimensions + # https://gis.stackexchange.com/questions/470207/how-to-write-crs-info-to-netcdf-in-a-way-qgis-can-read-python-xarray + + return ds + + open_dataset_parameters = ["filename", "drop_variables"] + + description = "Use Clawpack fixed grid output files in Xarray" + url = "https://www.clawpack.org/fgout.html" + + +class FGMaxBackend(BackendEntrypoint): + "Xarray Backend for Clawpack fgmax grid format." + + def open_dataset( + self, + filename, + epsg=None, + drop_variables=None, + clip=False, # if True, clip the entire array to the extent where arrival_time is defined. + ): + + if drop_variables is None: + drop_variables = [] + + # expectation is for standard clawpack organization, + # e.g., output and 'fgmax_grids.data' in _output + fgno = int(os.path.basename(filename).split(".")[0][-4:]) + outdir = os.path.dirname(filename) + data_file = os.path.join(outdir, "fgmax_grids.data") + + fg = fgmax_tools.FGmaxGrid() + fg.read_fgmax_grids_data(fgno=fgno, data_file=data_file) + if fg.point_style != 2: + raise ValueError("FGMaxBackend only works with fg.point_style=2") + + fg.read_output(outdir=outdir) + + # Construct the x and y coordinates + # Both come in ascending, therefore flip y so that it is ordered as expected. + x = fg.x + y = np.flipud(fg.y) + + # Construct the data_vars array. To organize the + # data in the way expected by netcdf standards, need + # to both transpose and flipud the array. + + data_vars = {} + + data_vars["arrival_time"] = _prepare_var( + fg.arrival_time.data, + fg.arrival_time.mask, + [ + "y", + "x", + ], + "seconds", + nodata, + "Wave arrival time", + ) + + data_vars["h_max"] = _prepare_var( + fg.h.data, + fg.h.mask, + [ + "y", + "x", + ], + "meters", + nodata, + "Maximum water depth", + ) + + data_vars["eta_max"] = _prepare_var( + fg.h.data + fg.B.data, + fg.h.mask, + [ + "y", + "x", + ], + "meters", + nodata, + "Maximum water surface elevation", + ) + + data_vars["h_max_time"] = _prepare_var( + fg.h_time.data, + fg.h_time.mask, + [ + "y", + "x", + ], + "seconds", + nodata, + "Time of maximum water depth", + ) + + data_vars["B"] = _prepare_var( + fg.B, + None, + [ + "y", + "x", + ], + "meters", + nodata, + "Basal topography at the first time fgmax first monitored maximum amr level", + ) + + data_vars["level"] = _prepare_var( + fg.level, + None, + [ + "y", + "x", + ], + "(no units)", + -1, + "Maximum amr level", + ) + + if hasattr(fg, "s"): + if hasattr(fg.s, "data"): + data_vars["s_max"] = _prepare_var( + fg.s.data, + fg.s.mask, + [ + "y", + "x", + ], + "meters per second", + nodata, + "Maximum velocity", + ) + data_vars["s_max_time"] = _prepare_var( + fg.s_time.data, + fg.s_time.mask, + [ + "y", + "x", + ], + "seconds", + nodata, + "Time of maximum velocity", + ) + if hasattr(fg, "hs"): + if hasattr(fg.hs, "data"): + data_vars["hs_max"] = _prepare_var( + fg.hs.data, + fg.hs.mask, + [ + "y", + "x", + ], + "meters squared per second", + nodata, + "Maximum momentum", + ) + data_vars["hs_max_time"] = _prepare_var( + fg.hs_time.data, + fg.hs_time.mask, + [ + "y", + "x", + ], + "seconds", + nodata, + "Time of maximum momentum", + ) + + data_vars["hss_max"] = _prepare_var( + fg.hss.data, + fg.hss.mask, + [ + "y", + "x", + ], + "meters cubed per second squared", + nodata, + "Maximum momentum flux", + ) + data_vars["hss_max_time"] = _prepare_var( + fg.hss_time.data, + fg.hss_time.mask, + [ + "y", + "x", + ], + "seconds", + nodata, + "Time of maximum momentum flux", + ) + + data_vars["h_min"] = _prepare_var( + fg.hmin.data, + fg.hmin.mask, + [ + "y", + "x", + ], + "meters", + nodata, + "Minimum depth", + ) + data_vars["h_min_time"] = _prepare_var( + fg.hmin_time.data, + fg.hmin_time.mask, + [ + "y", + "x", + ], + "seconds", + nodata, + "Time of minimum depth", + ) + + # drop requested variables + for var in drop_variables: + if var in data_vars.keys(): + del data_vars[var] + + # Construct the values from + ds_attrs = {"description": "D-Claw model output"} + + ds = xr.Dataset( + data_vars=data_vars, + coords=dict( + x=(["x"], x, {"units": "meters"}), + y=(["y"], y, {"units": "meters"}), + ), + attrs=ds_attrs, + ) + + if epsg is not None: + + ds.rio.write_crs( + epsg, + inplace=True, + ).rio.set_spatial_dims( + x_dim="x", + y_dim="y", + inplace=True, + ).rio.write_coordinate_system(inplace=True) + # https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Spatial-dimensions + # https://gis.stackexchange.com/questions/470207/how-to-write-crs-info-to-netcdf-in-a-way-qgis-can-read-python-xarray + + # clip + if clip: + clip_data = fg.arrival_time.data >= 0 + clip_data = clip_data.T + clip_data = np.flipud(clip_data) + ds = _clip(ds, clip_data) + + return ds + + open_dataset_parameters = ["filename", "drop_variables", "clip"] + + description = "Use Clawpack fix grid monitoring files in Xarray" + url = "https://www.clawpack.org/fgmax.html#fgmax" + + +def _clip(ds, clip_data): + # clip DS to the extent where clip_data is true. + # useful for when there are large areas where nothing happened. + + nrow, ncol = clip_data.shape + + valid_col = np.sum(clip_data, axis=0) > 0 + valid_row = np.sum(clip_data, axis=1) > 0 + + sel_rows = np.nonzero(valid_row)[0] + sel_cols = np.nonzero(valid_col)[0] + + first_row = sel_rows[0] + last_row = sel_rows[-1] + first_col = sel_cols[0] + last_col = sel_cols[-1] + + ds = ds.isel(y=np.arange(first_row, last_row), x=np.arange(first_col, last_col)) + return ds