Merge branch 'master' into add-track-plotting

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.
+

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

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
 
 
     #-----------------------------------------

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

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

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

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")

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

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