[GRASS-dev] Fwd: [GRASS GIS] #3246: r.in.srtm: patch to support SRTM water bodies import (SWBD)

Luca Delucchi lucadeluge at gmail.com
Fri Feb 10 08:29:33 PST 2017


I'm not able to upload the new patch because there is an error [0].

I'm attaching the new patch here


[0] https://trac.osgeo.org/osgeo/ticket/1875

---------- Forwarded message ----------
From: GRASS GIS <trac at osgeo.org>
Date: 10 February 2017 at 17:17
Subject: Re: [GRASS GIS] #3246: r.in.srtm: patch to support SRTM water
bodies import (SWBD)
To:


#3246: r.in.srtm: patch to support SRTM water bodies import (SWBD)
--------------------------+---------------------------------
  Reporter:  stjo         |      Owner:  grass-dev@…
      Type:  enhancement  |     Status:  new
  Priority:  normal       |  Milestone:  7.2.1
 Component:  Python       |    Version:  svn-releasebranch72
Resolution:               |   Keywords:  r.in.srtm
       CPU:  Unspecified  |   Platform:  All
--------------------------+---------------------------------

Comment (by lucadelu):

 I updated the patch, but I'm not able to test it because there are no raw
 files in https://dds.cr.usgs.gov/srtm/version2_1/SWBD/. In the zip files
 there are only shapefile.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3246#comment:1>
GRASS GIS <https://grass.osgeo.org>



-- 
ciao
Luca

www.lucadelu.org
-------------- next part --------------
Index: r.in.srtm.html
===================================================================
--- r.in.srtm.html	(revision 70529)
+++ r.in.srtm.html	(working copy)
@@ -1,41 +1,73 @@
 <h2>DESCRIPTION</h2>
 
-<em>r.in.srtm</em> imports SRTM hgt files into GRASS.
+<em>r.in.srtm</em> imports SRTM hgt files or SRTM SWBD (SRTM Water Body Data) raw
+files into GRASS GIS.
 
-SRTM Version 1 and improved Version 2 data sets can be downloaded from 
-NASA at this site:<br>
-<a href="http://dds.cr.usgs.gov/srtm/">http://dds.cr.usgs.gov/srtm/</a>
+SRTM Version 1 and improved Version 2 data sets can be downloaded from
+NASA at <a href="http://dds.cr.usgs.gov/srtm/">this site</a>.<br>
 
 <p>
-Gap-filled SRTM Version 3 data can be downloaded from USGS at this site:<br>
-<a href="http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL3.003/2000.02.11/">http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL3.003/2000.02.11/</a>
+Gap-filled SRTM Version 3 (1 arc seconds at 30m) data can be downloaded from USGS
+at <a href="http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL3.003/2000.02.11/">this site</a>.
 
+<p>
+SRTM SWBD (SRTM Water Body Data; 1 arc seconds at 30m) can be downloaded from USGS at
+<a href="http://e4ftl01.cr.usgs.gov/SRTM/SRTMSWBD.003/2000.02.11/">this site</a>.
+
+<p>
+The default output is the input name, without extension. It is also possible to
+import the zip files directly.
+
 <h2>NOTES</h2>
 
 SRTM tiles are of 1 degree by 1 degree size. The SRTM filename contains the
-coordinates which refer to the <b>center</b> of the lower left pixel (e.g., N51E010: 
+coordinates which refer to the <b>center</b> of the lower left pixel (e.g., N51E010:
 lower left cell center at 10E, 51N). To identify a tile name, a grid can be easily
 visualized in the GRASS monitor:
 
 <div class="code"><pre>
+# in latitude-longitude location
+d.mon wx0
 d.grid size=1
 </pre></div>
 
-To import TOPEX/SRTM30 PLUS data, use <em>r.in.bin</em>.
+Hint: To import TOPEX/SRTM30 PLUS data, use <em>r.in.bin</em>.
 
+<h2>EXAMPLES</h2>
+
+<h3>Import of SRTM</h3>
+Import of SRTM V3 (1 arc seconds at 30m):
+<div class="code"><pre>
+r.in.srtm -1 input=N00E006.hgt output=N00E006
+</pre></div>
+
+Import of SRTM V3 zip file (1 arc seconds at 30m):
+<div class="code"><pre>
+r.in.srtm -1 input=N00E006.SRTMGL1.hgt.zip output=N00E006
+</pre></div>
+
+<h3>Import of SRTM water bodies</h3>
+
+Import of SRTM SWBD (1 arc seconds at 30m):
+<div class="code"><pre>
+# using predefined output
+r.in.srtm -1 -w input=N00E006.raw
+</pre></div>
+
 <h2>SEE ALSO</h2>
 
 <em>
-<a href="r.in.bin.html">r.in.bin</a>
+<a href="r.import.html">r.import</a>,
+<a href="r.in.bin.html">r.in.bin</a>,
+<a href="r.in.gdal.html">r.in.gdal</a>
 </em>
 <p>The <a href="http://www2.jpl.nasa.gov/srtm/">Shuttle Radar Topography Mission</a>
 homepage at NASA's JPL.
-<br>
-The <a href="http://pub7.bravenet.com/forum/537683448/">SRTM Web Forum</a>
 
 <h2>REFERENCES</h2>
 
-M. Neteler, 2005. <a href="http://grass.osgeo.org/newsletter/GRASSNews_vol3.pdf">SRTM and VMAP0 data in OGR and GRASS.</a> <i><a href="http://grass.osgeo.org/newsletter/">GRASS Newsletter</a></i>, Vol.3, pp. 2-6, June 2005. ISSN 1614-8746.
+M. Neteler, 2005. <a href="http://grass.osgeo.org/newsletter/GRASSNews_vol3.pdf">SRTM and VMAP0 data in OGR and GRASS.</a>
+<i><a href="http://grass.osgeo.org/newsletter/">GRASS Newsletter</a></i>, Vol.3, pp. 2-6, June 2005. ISSN 1614-8746.
 
 
 <h2>AUTHORS</h2>
@@ -42,6 +74,7 @@
 
 Markus Neteler<br>
 Improved by W. Kyngesburye and H. Bowman<br>
-Update for SRTM V3 by Markus Metz
+Update for SRTM V3 by Markus Metz<br>
+Update for SRTM SWBD by Jonas Strobel
 
 <p><i>Last changed: $Date$</i>
Index: r.in.srtm.py
===================================================================
--- r.in.srtm.py	(revision 70529)
+++ r.in.srtm.py	(working copy)
@@ -2,13 +2,14 @@
 #
 ############################################################################
 #
-# MODULE:    r_in_aster.py
+# MODULE:    r.in.srtm.py
 # AUTHOR(S): Markus Neteler 11/2003 neteler AT itc it
 #            Hamish Bowman
 #            Glynn Clements
-# PURPOSE:   import of SRTM hgt files into GRASS
+#            Jonas Strobel
+# PURPOSE:   import of SRTM hgt files or SRTM SWBD raw files into GRASS
 #
-# COPYRIGHT:	(C) 2004, 2006 by the GRASS Development Team
+# COPYRIGHT:    (C) 2004-2017 by the GRASS Development Team
 #
 #		This program is free software under the GNU General Public
 #		License (>=v2). Read the file COPYING that comes with GRASS
@@ -22,6 +23,7 @@
 # April 2006: links updated from ftp://e0dps01u.ecs.nasa.gov/srtm/
 #             to current links below
 # October 2008: Converted to Python by Glynn Clements
+# December 2016: Added flag to read SRTM SWBD (SRTM Water Body Dataset) raw binary files
 #########################
 # Derived from:
 # ftp://e0srp01u.ecs.nasa.gov/srtm/version1/Documentation/Notes_for_ARCInfo_users.txt
@@ -47,14 +49,22 @@
 # http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/intpt.html
 #
 #########################
+# Test cases
+# - SRTM HGT 1 arc second
+# - SRTM HGT 3 arc second
+# - SRTM HGT.zip 1 arc second
+# - SRTM HGT.zip 3 arc second
+# - SRTM SWBD 1 arc second
+# - SRTM SWBD.zip 1 arc second
+#########################
 
 #%Module
-#% description: Imports SRTM HGT files into raster map.
+#% description: Imports SRTM HGT or RAW files into raster map.
 #% keyword: raster
 #% keyword: import
 #%End
 #%option G_OPT_F_INPUT
-#% description: Name of SRTM input tile (file without .hgt.zip extension)
+#% description: Name of SRTM HGT / SRTM RAW input tile
 #%end
 #%option G_OPT_R_OUTPUT
 #% description: Name for output raster map (default: input tile)
@@ -64,6 +74,10 @@
 #% key: 1
 #% description: Input is a 1-arcsec tile (default: 3-arcsec)
 #%end
+#%flag
+#% key: w
+#% description: Import SRTM SWBD (SRTM Water Body Data)
+#%end
 
 tmpl1sec = """BYTEORDER M
 LAYOUT BIL
@@ -117,7 +131,7 @@
 def cleanup():
     if not in_temp:
         return
-    for ext in ['.bil', '.hdr', '.prj', '.hgt.zip']:
+    for ext in ['.bil', '.hdr', '.prj', '.hgt.zip', '.raw', '.raw.zip']:
         grass.try_remove(tile + ext)
     os.chdir('..')
     grass.try_rmdir(tmpdir)
@@ -131,6 +145,7 @@
     input = options['input']
     output = options['output']
     one = flags['1']
+    water = flags['w']
 
     # are we in LatLong location?
     s = grass.read_command("g.proj", flags='j')
@@ -140,17 +155,23 @@
 
     # use these from now on:
     infile = input
-    while infile[-4:].lower() in ['.hgt', '.zip']:
+    while infile[-4:].lower() in ['.hgt', '.zip', '.raw']:
         infile = infile[:-4]
     (fdir, tile) = os.path.split(infile)
 
     if not output:
         tileout = tile
+        grass.debug("No output set... using name: " + tileout)
     else:
         tileout = output
 
-    zipfile = infile + ".hgt.zip"
-    hgtfile = os.path.join(fdir, tile[:7] + ".hgt")
+    if not water:
+        zipfile = infile + ".hgt.zip"
+        hgtfile = os.path.join(fdir, tile[:7] + ".hgt")
+    else:
+        zipfile = infile + ".raw.zip"
+        hgtfile = os.path.join(fdir, tile[:7] + ".raw")
+
     if os.path.isfile(zipfile):
         # check if we have unzip
         if not grass.find_program('unzip'):
@@ -162,11 +183,16 @@
         tenv['UNZIP'] = '-qq'
         if grass.call(['unzip', '-t', zipfile], env=tenv) != 0:
             grass.fatal(_("'%s' does not appear to be a valid zip file.") % zipfile)
+        is_zip = True
+    elif not water:
+        os.path.isfile(hgtfile)
+        # try and see if it's already unzipped
+        is_zip = False
 
-        is_zip = True
     elif os.path.isfile(hgtfile):
         # try and see if it's already unzipped
         is_zip = False
+
     else:
         grass.fatal(_("File '%s' or '%s' not found") % (zipfile, hgtfile))
 
@@ -176,26 +202,39 @@
     os.mkdir(tmpdir)
 
     if is_zip:
-        shutil.copyfile(zipfile, os.path.join(tmpdir, tile + ".hgt.zip"))
+        if not water:
+            shutil.copyfile(zipfile, os.path.join(tmpdir, tile + ".hgt.zip"))
+        else:
+            shutil.copyfile(zipfile, os.path.join(tmpdir, tile + ".raw.zip"))
+
     else:
-        shutil.copyfile(hgtfile, os.path.join(tmpdir, tile + ".hgt"))
+        if not water:
+            shutil.copyfile(hgtfile, os.path.join(tmpdir, tile + ".hgt"))
+        else:
+            shutil.copyfile(hgtfile, os.path.join(tmpdir, tile + ".raw"))
 
     # change to temporary directory
     os.chdir(tmpdir)
     in_temp = True
 
-    zipfile = tile + ".hgt.zip"
-    hgtfile = tile[:7] + ".hgt"
+    if not water:
+        zipfile = tile + ".hgt.zip"
+        hgtfile = tile[:7] + ".hgt"
+    else:
+        zipfile = tile + ".raw.zip"
+        hgtfile = tile[:7] + ".raw"
+
     bilfile = tile + ".bil"
 
     if is_zip:
         # unzip & rename data file:
-        grass.message(_("Extracting '%s'...") % infile)
+        grass.message(_("Extracting '%s'...") % zipfile)
         if grass.call(['unzip', zipfile], env=tenv) != 0:
             grass.fatal(_("Unable to unzip file."))
 
-    grass.message(_("Converting input file to BIL..."))
-    os.rename(hgtfile, bilfile)
+    if not water:
+        grass.message(_("Converting input file to BIL..."))
+        os.rename(hgtfile, bilfile)
 
     north = tile[0]
     ll_latitude = int(tile[1:3])
@@ -210,17 +249,27 @@
     if east == "W":
         ll_longitude *= -1
 
-    # Calculate Upper Left from Lower Left
-    ulxmap = "%.1f" % ll_longitude
-    # SRTM90 tile size is 1 deg:
-    ulymap = "%.1f" % (ll_latitude + 1)
+    if water:
+        # Calculate Upper Left from Lower Left
+        ulxmap = "%.1f" % (ll_longitude + 1)
 
-    if not one:
-        tmpl = tmpl3sec
+        # SRTM90 tile size is 1 deg:
+        ulymap = "%.1f" % (ll_latitude + 1)
+
     else:
+        # Calculate Upper Left from Lower Left
+        ulxmap = "%.1f" % ll_longitude
+
+        # SRTM90 tile size is 1 deg:
+        ulymap = "%.1f" % (ll_latitude + 1)
+
+    if one or water:
         grass.message(_("Attempting to import 1-arcsec data."))
         tmpl = tmpl1sec
 
+    else:
+        tmpl = tmpl3sec
+
     header = tmpl % (ulxmap, ulymap)
     hdrfile = tile + '.hdr'
     outf = file(hdrfile, 'w')
@@ -233,20 +282,38 @@
     outf.write(proj)
     outf.close()
 
-    try:
-        grass.run_command('r.in.gdal', input=bilfile, out=tileout)
-    except:
-        grass.fatal(_("Unable to import data"))
+    if not water:
+        try:
+            grass.run_command('r.in.gdal', input=bilfile, output=tileout)
+        except:
+            grass.fatal("Unable to import data")
 
+    else:
+        # If water, these operations are required
+        swbd_res = 0.000277777777777778  # 0:00:01
+
+        n = float(ulymap) + (0.5 * swbd_res)
+        s = float(ll_latitude) - (0.5 * swbd_res)
+        e = float(ulxmap) + (0.5 * swbd_res)
+        w = int(ll_longitude) - (0.5 * swbd_res)
+
+        try:
+            grass.run_command('r.in.bin', input=hgtfile, output=tileout,
+                              bytes=1, north=n, south=s, east=e, west=w,
+                              rows=3601, cols=3601)
+        except:
+            grass.fatal(_("Unable to import data"))
+
     # nice color table
     grass.run_command('r.colors', map=tileout, color='srtm')
 
     # write cmd history:
     grass.raster_history(tileout)
-
     grass.message(_("Done: generated map ") + tileout)
-    grass.message(_("(Note: Holes in the data can be closed with 'r.fillnulls' using splines)"))
 
+    if not water:
+        grass.message(_("(Note: Holes in the data can be closed with 'r.fillnulls' using splines)"))
+
 if __name__ == "__main__":
     options, flags = grass.parser()
     atexit.register(cleanup)


More information about the grass-dev mailing list