Merge pull request #210 from rjleveque/topotool_updates

WIP: Topotool updates

examples/tsunami/bowl-radial/maketopo.py

@@ -3,7 +3,7 @@
 Module to create topo and qinit data files for this example.
 """
 
-from clawpack.geoclaw import topotools
+from clawpack.geoclaw.topotools import Topography
 from numpy import *
 
 def maketopo():
@@ -17,7 +17,11 @@ def maketopo():
     yupper = 100.e0
     ylower = -100.e0
     outfile= "bowl.topotype2"     
-    topotools.topo2writer(outfile,topo,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+
+    topography = Topography(topo_func=topo)
+    topography.x = linspace(xlower,xupper,nxpoints)
+    topography.y = linspace(ylower,yupper,nypoints)
+    topography.write(outfile, topo_type=2, Z_format="%22.15e")
 
 def makeqinit():
     """
@@ -30,7 +34,11 @@ def makeqinit():
     yupper = 50.e0
     ylower = -50.e0
     outfile= "hump.xyz"     
-    topotools.topo1writer(outfile,qinit,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+
+    topography = Topography(topo_func=qinit)
+    topography.x = linspace(xlower,xupper,nxpoints)
+    topography.y = linspace(ylower,yupper,nypoints)
+    topography.write(outfile, topo_type=1)
 
 def topo(x,y):
     """

examples/tsunami/bowl-slosh/maketopo.py

@@ -3,7 +3,7 @@
 Module to create topo and qinit data files for this example.
 """
 
-from clawpack.geoclaw.topotools import topo1writer, topo2writer
+from clawpack.geoclaw.topotools import Topography
 from numpy import *
 
 #from pyclaw.data import Data
@@ -26,7 +26,12 @@ def maketopo():
     xlower = -2.e0
     ylower = -2.e0
     outfile= "bowl.topotype2"
-    topo2writer(outfile,topo,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+
+    topography = Topography(topo_func=topo)
+    topography.x = linspace(xlower,xupper,nxpoints)
+    topography.y = linspace(ylower,yupper,nypoints)
+    topography.write(outfile, topo_type=2, Z_format="%22.15e")
+
 
 def topo(x,y):
     """

src/python/geoclaw/etopotools.py

@@ -59,7 +59,7 @@ def etopo1_download(xlimits, ylimits, dx=0.0166666666667, dy=None, \
     y1,y2 = ylimits
 
     if file_name is None:
-        file_name = 'etopo1_%i_%i_%i_%i_%imin.tt3' \
+        file_name = 'etopo1_%i_%i_%i_%i_%imin.asc' \
             % (int(round(x1)),int(round(x2)),int(round(y1)),int(round(y2)),\
               int(round(60*dx)))
 

src/python/geoclaw/topotools.py

@@ -788,8 +788,8 @@ def read_header(self):
 
         return num_cells
 
-    def write(self, path, no_data_value=None, topo_type=None, masked=True, 
-                header_style='geoclaw'):
+    def write(self, path, topo_type=None, no_data_value=None, masked=True, 
+                header_style='geoclaw', Z_format="%15.7e"):
         r"""Write out a topography file to path of type *topo_type*.
 
         Writes out a topography file of topo type specified with *topo_type* or
@@ -799,13 +799,25 @@ def write(self, path, no_data_value=None, topo_type=None, masked=True,
 
         :Input:
          - *path* (str)  - file to write
-         - *no_data_value* - values used to indicate missing data
          - *topo_type* (int) - GeoClaw format topo_type 
+           **Note:** this is second positional argument, agreeing with
+           the read function in this class.  It was the third argument in
+           GeoClaw version 5.3.1 and earlier.  
+         - *no_data_value* - values used to indicate missing data
          - *masked* (bool) - unused??
          - *header_style* (str) - indicates format of header lines
              'geoclaw' or 'default'  ==> write value then label 
              'arcgis' or 'asc' ==> write label then value  
                         (needed for .asc files in ArcGIS)
+         - *Z_format* (str) - string format to use for Z values
+           The default format "%15.7e" gives at least millimeter precision
+           for topography with abs(Z) < 10000 and results in
+           smaller files than the previous default of "%22.15e" used in
+           GeoClaw version 5.3.1 and earlier.  A shorter format can be used
+           if the user knows there are fewer significant digits, e.g.
+           etopo1 data is integers and so has a resolution of 1 meter.
+           In this case a cropped or coarsened version might be written
+           with `Z_format = "%7i"`, for example.
 
         """
 
@@ -883,23 +895,25 @@ def write(self, path, no_data_value=None, topo_type=None, masked=True,
 
                 # Write out topography data
                 if topo_type == 2:
+                    Z_format = Z_format + "\n"
                     if masked_Z:
                         Z_filled = numpy.flipud(self.Z.filled())
                     else:
                         Z_filled = numpy.flipud(self.Z)
                     for i in xrange(self.Z.shape[0]):
                         for j in xrange(self.Z.shape[1]):
-                            outfile.write("%22.15e\n" % Z_filled[i,j])
+                            outfile.write(Z_format % Z_filled[i,j])
                     if masked_Z:
                         del Z_filled
                 elif topo_type == 3:
+                    Z_format = Z_format + " "
                     if masked_Z:
                         Z_flipped = numpy.flipud(self.Z.filled())
                     else:
                         Z_flipped = numpy.flipud(self.Z)
                     for i in xrange(self.Z.shape[0]):
                         for j in xrange(self.Z.shape[1]):
-                            outfile.write("%22.15e   " % (Z_flipped[i,j]))
+                            outfile.write(Z_format % (Z_flipped[i,j]))
                         outfile.write("\n")
                     if masked_Z:
                         del Z_flipped
@@ -1414,7 +1428,7 @@ def smooth_data(self, indices, r=1):
                 self.Z[index[0], index[1]] = summation / num_points
 
 
-    def crop(self, filter_region):
+    def crop(self, filter_region=None, coarsen=1):
         r"""Crop region to *filter_region*
 
         Create a new Topography object that is identical to this one but cropped
@@ -1429,6 +1443,10 @@ def crop(self, filter_region):
         if self.unstructured:
             raise NotImplemented("*** Cannot currently crop unstructured topo")
 
+        if filter_region is None:
+            # only want to coarsen, so this is entire region:
+            filter_region = [self.x[0],self.x[-1],self.y[0],self.y[-1]]
+
         # Find indices of region
         region_index = [None, None, None, None]
         region_index[0] = (self.x >= filter_region[0]).nonzero()[0][0]
@@ -1437,17 +1455,17 @@ def crop(self, filter_region):
         region_index[3] = (self.y <= filter_region[3]).nonzero()[0][-1] + 1
         newtopo = Topography()
 
-        newtopo._x = self._x[region_index[0]:region_index[1]]
-        newtopo._y = self._y[region_index[2]:region_index[3]]
+        newtopo._x = self._x[region_index[0]:region_index[1]:coarsen]
+        newtopo._y = self._y[region_index[2]:region_index[3]:coarsen]
 
         # Force regeneration of 2d coordinate arrays and extent if needed
         newtopo._X = None
         newtopo._Y = None
         newtopo._extent = None
 
         # Modify Z array as well
-        newtopo._Z = self._Z[region_index[2]:region_index[3],
-                          region_index[0]:region_index[1]]
+        newtopo._Z = self._Z[region_index[2]:region_index[3]:coarsen,
+                          region_index[0]:region_index[1]:coarsen]
 
         newtopo.unstructured = self.unstructured
         newtopo.topo_type = self.topo_type

tests/bowl_slosh/maketopo.py

@@ -3,7 +3,8 @@
 Module to create topo and qinit data files for this example.
 """
 
-from clawpack.geoclaw.topotools import topo1writer, topo2writer
+from clawpack.geoclaw.topotools import Topography
+
 from numpy import *
 
 #from pyclaw.data import Data
@@ -26,7 +27,10 @@ def maketopo():
     xlower = -2.e0
     ylower = -2.e0
     outfile= "bowl.topotype2"
-    topo2writer(outfile,topo,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+    topography = Topography(topo_func=topo)
+    topography.x = linspace(xlower,xupper,nxpoints)
+    topography.y = linspace(ylower,yupper,nypoints)
+    topography.write(outfile, topo_type=2, Z_format="%22.15e")
 
 def topo(x,y):
     """

tests/bowl_slosh/regression_tests.py

@@ -33,7 +33,8 @@ def setUp(self):
         topo.topo_type = 2
         topo.x = numpy.linspace(-2.0, 2.0, 200)
         topo.y = numpy.linspace(-2.0, 2.0, 200)
-        topo.write(os.path.join(self.temp_path, "bowl.topotype2"))
+        topo.write(os.path.join(self.temp_path, "bowl.topotype2"), \
+                topo_type=2, Z_format="%22.15e")
 
         from make_fgmax_grid import make_fgmax_grid1
         make_fgmax_grid1(self.temp_path)

tests/dtopo1/maketopo.py

@@ -3,7 +3,8 @@
 Module to create topo and qinit data files for this example.
 """
 
-from clawpack.geoclaw.topotools import topo1writer, topo2writer
+from clawpack.geoclaw.topotools import Topography
+
 from clawpack.geoclaw import dtopotools
 import numpy
 import matplotlib.pyplot as plt
@@ -21,7 +22,10 @@ def maketopo():
     ylower = -10e0
     yupper= 10e0
     outfile= "topo1.topotype2"
-    topo2writer(outfile,topo,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+    topography = Topography(topo_func=topo)
+    topography.x = numpy.linspace(xlower,xupper,nxpoints)
+    topography.y = numpy.linspace(ylower,yupper,nypoints)
+    topography.write(outfile, topo_type=2, Z_format="%22.15e")
     dx = (xupper-xlower)/(nxpoints-1)
     dy = (yupper-ylower)/(nypoints-1)
     print "==> topo1 has dx = %g, dy = %g" % (dx,dy)
@@ -37,7 +41,10 @@ def maketopo2():
     ylower = -0.1
     yupper= 0.4
     outfile= "topo2.topotype2"
-    topo2writer(outfile,topo2,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+    topography = Topography(topo_func=topo2)
+    topography.x = numpy.linspace(xlower,xupper,nxpoints)
+    topography.y = numpy.linspace(ylower,yupper,nypoints)
+    topography.write(outfile, topo_type=2, Z_format="%22.15e")
     dx = (xupper-xlower)/(nxpoints-1)
     dy = (yupper-ylower)/(nypoints-1)
     print "==> topo2 has dx = %g, dy = %g" % (dx,dy)