[GRASS-SVN] r50079 - grass-addons/grass7/raster/r.in.srtm.region
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jan 7 06:30:50 EST 2012
Author: mmetz
Date: 2012-01-07 03:30:50 -0800 (Sat, 07 Jan 2012)
New Revision: 50079
Added:
grass-addons/grass7/raster/r.in.srtm.region/Makefile
grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.html
grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.py
Log:
import new module r.in.srtm.region
Added: grass-addons/grass7/raster/r.in.srtm.region/Makefile
===================================================================
--- grass-addons/grass7/raster/r.in.srtm.region/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.in.srtm.region/Makefile 2012-01-07 11:30:50 UTC (rev 50079)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.in.srtm.region
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.html
===================================================================
--- grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.html (rev 0)
+++ grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.html 2012-01-07 11:30:50 UTC (rev 50079)
@@ -0,0 +1,36 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.in.srtm.region</em> imports all SRTM V2.1 tiles covering the
+current region or region extents given with <b>region</b> into GRASS,
+patches the tiles together and optionally interpolates holes.
+<p>
+<em>r.in.srtm.region</em> downloads SRTM tiles from:<br>
+<a href="http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/">http://dds.cr.usgs.gov/srtm/</a>
+<br>
+or optionally uses a local folder with previously downloaded files if
+the <b>local</b> option is given.
+
+<h2>NOTES</h2>
+
+<em>r.in.srtm.region</em> imports only SRTM tiles with 3 arcsec resolution.
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.in.srtm.html">r.in.srtm</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.
+
+
+<h2>AUTHORS</h2>
+
+Markus Metz<br>
+
+<p><i>Last changed: $Date: 2011-11-18 08:53:00 +0100 (Fri, 18 Nov 2011) $</i>
Added: grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.py
===================================================================
--- grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.py (rev 0)
+++ grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.py 2012-01-07 11:30:50 UTC (rev 50079)
@@ -0,0 +1,347 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: r.in.srtm.region
+#
+# AUTHOR(S): Markus Metz
+#
+# PURPOSE: Create a DEM from 3 arcsec SRTM v2.1 tiles
+#
+# COPYRIGHT: (C) 2011 GRASS development team
+#
+# This program is free software under the GNU General
+# Public License (>=v2). Read the file COPYING that
+# comes with GRASS for details.
+#
+#############################################################################
+
+#%module
+#% description: Creates a DEM from 3 arcsec SRTM v2.1 tiles.
+#% keywords: raster
+#% keywords: import
+#%end
+#%option
+#% key: url
+#% description: base url to fetch SRTM v2.1 tiles
+#% answer: http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/
+#% required: no
+#%end
+#%option G_OPT_M_DIR
+#% key: local
+#% label: local folder with SRTM v2.1 tiles
+#% description: use local folder instead of url to retrieve SRTM tiles
+#% required: no
+#%end
+#%option
+#% key: region
+#% type: double
+#% label: Import subregion only (default is current region)
+#% description: Format: xmin,ymin,xmax,ymax - usually W,S,E,N
+#% key_desc: xmin,ymin,xmax,ymax
+#% multiple: yes
+#% required: no
+#%end
+#%option G_OPT_R_OUTPUT
+#% description: Name for output raster map
+#% required: yes
+#%end
+#%option
+#% key: memory
+#% type: integer
+#% description: Memory in MB for interpolation
+#% answer: 300
+#% required: no
+#%end
+#%flag
+#% key: n
+#% description: Fill null cells
+#%end
+
+
+proj = ''.join([
+ 'GEOGCS[',
+ '"wgs84",',
+ 'DATUM["WGS_1984",SPHEROID["wgs84",6378137,298.257223563],TOWGS84[0.000000,0.000000,0.000000]],',
+ 'PRIMEM["Greenwich",0],',
+ 'UNIT["degree",0.0174532925199433]',
+ ']'])
+
+import sys
+import os
+import shutil
+import atexit
+import urllib
+import urllib2
+import time
+
+import grass.script as grass
+
+def import_local_tile(tile, local, pid):
+ output = tile + '.r.in.srtm2.tmp.' + str(pid)
+ local_tile = str(tile) + '.hgt.zip'
+
+ path = os.path.join(local, local_tile)
+ if os.path.exists(path):
+ path = os.path.join(local, tile)
+ grass.run_command('r.in.srtm', input = path, output = output, quiet = True)
+ return 1
+
+ # SRTM subdirs: Africa, Australia, Eurasia, Islands, North_America, South_America
+ for srtmdir in ('Africa', 'Australia', 'Eurasia', 'Islands', 'North_America', 'South_America'):
+ path = os.path.join(local, srtmdir, local_tile)
+
+ if os.path.exists(path):
+ path = os.path.join(local, srtmdir, tile)
+ grass.run_command('r.in.srtm', input = path, output = output, quiet = True)
+ return 1
+
+ return 0
+
+def download_tile(tile, url, pid):
+ output = tile + '.r.in.srtm2.tmp.' + str(pid)
+ local_tile = str(tile) + '.hgt.zip'
+
+ urllib.urlcleanup()
+
+ # SRTM subdirs: Africa, Australia, Eurasia, Islands, North_America, South_America
+ for srtmdir in ('Africa', 'Australia', 'Eurasia', 'Islands', 'North_America', 'South_America'):
+ remote_tile = str(url) + str(srtmdir) + '/' + local_tile
+ goturl = 1
+
+ try:
+ f = urllib2.urlopen(remote_tile)
+ fo = open(local_tile, 'w+b')
+ fo.write(f.read())
+ fo.close
+ time.sleep(0.5)
+ # does not work:
+ #urllib.urlretrieve(remote_tile, local_tile, data = None)
+ except:
+ goturl = 0
+ pass
+
+ if goturl == 1:
+ return 1
+
+ return 0
+
+
+def cleanup():
+ if not in_temp:
+ return
+ os.chdir(currdir)
+ grass.run_command('g.region', region = tmpregionname)
+ grass.run_command('g.remove', region = tmpregionname, quiet = True)
+ #grass.try_rmdir(tmpdir)
+
+def main():
+ global tile, tmpdir, in_temp, currdir, tmpregionname
+
+ in_temp = False
+
+ url = options['url']
+ local = options['local']
+ output = options['output']
+ memory = options['memory']
+ fillnulls = flags['n']
+
+ if len(url) == 0 and len(local) == 0:
+ grass.fatal(_("Either 'url' or 'local' is needed"))
+
+ if len(local) == 0:
+ local = None
+
+ # are we in LatLong location?
+ s = grass.read_command("g.proj", flags='j')
+ kv = grass.parse_key_val(s)
+ if kv['+proj'] != 'longlat':
+ grass.fatal(_("This module only operates in LatLong locations"))
+
+ if fillnulls == 1 and memory <= 0:
+ grass.warning(_("Amount of memory to use for interpolation must be positive, setting to 300 MB"))
+ memory = '300'
+
+ # make a temporary directory
+ tmpdir = grass.tempfile()
+ grass.try_remove(tmpdir)
+ os.mkdir(tmpdir)
+ currdir = os.getcwd()
+ pid = os.getpid()
+
+ # change to temporary directory
+ os.chdir(tmpdir)
+ in_temp = True
+ if local is None:
+ local = tmpdir
+
+ # get extents
+ reg = grass.region()
+ tmpregionname = 'r_in_srtm2_tmp_region'
+ grass.run_command('g.region', save = tmpregionname)
+ if options['region'] is None or options['region'] == '':
+ north = reg['n']
+ south = reg['s']
+ east = reg['e']
+ west = reg['w']
+ else:
+ west, south, east, north = options['region'].split(',')
+ west = float(west)
+ south = float(south)
+ east = float(east)
+ north = float(north)
+
+ # adjust extents to cover SRTM tiles: 1 degree bounds
+ tmpint = int(north)
+ if tmpint < north:
+ north = tmpint + 1
+ else:
+ north = tmpint
+ tmpint = int(south)
+ if tmpint > south:
+ south = tmpint - 1
+ else:
+ south = tmpint
+ tmpint = int(east)
+ if tmpint < east:
+ east = tmpint + 1
+ else:
+ east = tmpint
+ tmpint = int(west)
+ if tmpint > west:
+ west = tmpint - 1
+ else:
+ west = tmpint
+
+ rows = abs(north - south)
+ cols = abs(east - west)
+ ntiles = rows * cols
+ grass.message(_("Importing %d SRTM tiles...") % ntiles, flag = 'i')
+ counter = 1
+
+ srtmtiles = ''
+ valid_tiles = 0
+ for ndeg in range(south, north):
+ for edeg in range(west, east):
+ grass.percent(counter, ntiles, 1)
+ counter += 1
+ if ndeg < 0:
+ tile = 'S'
+ else:
+ tile = 'N'
+ tile = tile + '%02d' % abs(ndeg)
+ if edeg < 0:
+ tile = tile + 'W'
+ else:
+ tile = tile + 'E'
+ tile = tile + '%03d' % abs(edeg)
+ grass.debug("Tile: %s" % tile, debug = 1)
+
+ if local != tmpdir:
+ gotit = import_local_tile(tile, local, pid)
+ else:
+ gotit = download_tile(tile, url, pid)
+ if gotit == 1:
+ gotit = import_local_tile(tile, tmpdir, pid)
+ if gotit == 1:
+ grass.verbose(_("Tile %s successfully imported") % tile)
+ valid_tiles += 1
+ else:
+ # create tile with zeros
+ # north
+ if ndeg < 0:
+ tmpn = '%02d:59:58.5S' % (abs(ndeg) - 2)
+ else:
+ tmpn = '%02d:00:01.5N' % (ndeg + 1)
+ # south
+ if ndeg <= 0:
+ tmps = '%02d:00:01.5S' % abs(ndeg)
+ else:
+ tmps = '%02d:59:58.5N' % (ndeg - 1)
+ # east
+ if edeg < 0:
+ tmpe = '%03d:59:58.5W' % (abs(edeg) - 2)
+ else:
+ tmpe = '%03d:00:01.5E' % (edeg + 1)
+ # west
+ if edeg <= 0:
+ tmpw = '%03d:00:01.5W' % abs(edeg)
+ else:
+ tmpw = '%03d:59:58.5E' % (edeg - 1)
+
+ grass.run_command('g.region', n = tmpn, s = tmps, e = tmpe, w = tmpw, res = '00:00:03')
+ grass.run_command('r.mapcalc', expression = "%s = 0" % (tile + '.r.in.srtm2.tmp.' + str(pid)), quiet = True)
+ grass.run_command('g.region', region = tmpregionname)
+
+
+ # g.mlist with fs = comma does not work ???
+ srtmtiles = grass.read_command('g.mlist',
+ type = 'rast',
+ pattern = '*.r.in.srtm2.tmp.%d' % pid,
+ fs = 'newline',
+ quiet = True)
+
+ srtmtiles = srtmtiles.splitlines()
+ srtmtiles = ','.join(srtmtiles)
+
+ if valid_tiles == 0:
+ grass.run_command('g.remove', rast = str(srtmtiles), quiet = True)
+ grass.warning(_("No tiles imported"))
+ if local != tmpdir:
+ grass.fatal(_("Please check if local folder <%s> is correct.") % local)
+ else:
+ grass.fatal(_("Please check internet connection and if url <%s> is correct.") % url)
+
+ grass.run_command('g.region', rast = str(srtmtiles));
+
+ if fillnulls == 0:
+ grass.run_command('r.patch', input = srtmtiles, output = output)
+ else:
+ ncells = grass.region()['cells']
+ if long(ncells) > 1000000000:
+ grass.message(_("%s cells to interpolate, this will take some time") % str(ncells), flag = 'i')
+ grass.run_command('r.patch', input = srtmtiles, output = output + '.holes')
+ mapstats = grass.parse_command('r.univar', map = output + '.holes', flags = 'g', quiet = True)
+ if mapstats['null_cells'] == '0':
+ grass.run_command('g.rename', rast = '%s,%s' % (output + '.holes', output), quiet = True)
+ else:
+ grass.run_command('r.resamp.bspline',
+ input = output + '.holes',
+ output = output + '.interp',
+ se = '0.0025', sn = '0.0025',
+ method = 'bilinear',
+ memory = memory,
+ flags = 'n')
+ grass.run_command('r.patch',
+ input = '%s,%s' % (output + '.holes',
+ output + '.interp'),
+ output = output + '.float',
+ flags = 'z')
+ grass.run_command('r.mapcalc', expression = '%s = round(%s)' % (output, output + '.float'))
+ grass.run_command('g.remove',
+ rast = '%s,%s,%s' % (output + '.holes', output + '.interp', output + '.float'),
+ quiet = True)
+
+ grass.run_command('g.remove', rast = str(srtmtiles), quiet = True)
+
+ # nice color table
+ grass.run_command('r.colors', map = output, color = 'srtm', quiet = True)
+
+ # write metadata:
+ tmphist = grass.tempfile()
+ f = open(tmphist, 'w+')
+ f.write(os.environ['CMDLINE'])
+ f.close()
+ grass.run_command('r.support', map = output,
+ loadhistory = tmphist,
+ description = 'generated by r.in.srtm.region',
+ source1 = 'SRTM V2.1',
+ source2 = (local if local != tmpdir else url))
+ grass.try_remove(tmphist)
+
+ grass.message(_("Done: generated map <%s>") % output)
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ main()
Property changes on: grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.py
___________________________________________________________________
Added: svn:executable
+
More information about the grass-commit
mailing list