r.import copied from addons to trunk

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@65833 15284696-431f-4ddb-bdfa-cd5b030d7da7

scripts/Makefile

@@ -34,6 +34,7 @@ SUBDIRS = \
 	r.colors.stddev \
 	r.fillnulls \
 	r.grow \
+    r.import \
 	r.in.aster \
 	r.in.srtm \
 	r.in.wms \

scripts/r.import/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.import
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

scripts/r.import/r.import.html

@@ -0,0 +1,56 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.import</em> imports selected bands from a GDAL raster datasouce 
+into the current location and mapset. If the projection of the input 
+does not match the projection of the location, the input is reprojected 
+into the current location. If the projection of the input does match 
+the projection of the location, the input is imported directly with <a 
+href="r.in.gdal.html">r.in.gdal</a>.
+
+<h4>Resolution</h4>
+<em>r.import</em> reports the estimated target resolution for each 
+input band. The estimated resolution will usually be some floating 
+point number, e.g. 271.301. In case option <b>resolution</b> is set to
+<em>estimated</em> (default), this floating point number will be used
+as target resolution. Since the target resolution should be typically the rounded 
+estimated resolution, e.g. 250 or 300 instead of 271.301, flag <b>-e</b>
+can be used first to obtain the estimate without importing the raster bands.
+Then the desired resolution is set with option <b>resolution_value</b> 
+and option <b>resolution</b>=<em>value</em>.
+For latlong locations, the resolution might be set to arc seconds, e.g. 1, 3, 7.5, 
+15, and 30 arc seconds are commonly used resolutions.
+
+<h4>Resampling methods</h4>
+<b>nearest</b> is the simplest method and the only possible method for 
+categorical data.
+<p>
+<b>bilinear</b> does linear interpolation and provides smoother output 
+than <b>nearest</b>. <b>bilinear</b> is recommended when reprojecting a 
+DEM for hydrological analysis or for surfaces where overshoots must be 
+avoided, e.g. precipitation should not become negative. 
+<p> 
+<b>bicubic</b> produces smoother output than <b>bilinear</b>, at 
+the cost of overshoots. 
+<p>
+<b>lanczos</b> produces the smoothest output of all methods and 
+preserves contrast best. <b>lanczos</b> is recommended for imagery. 
+Both <b>bicubic</b> and <b>lanczos</b> preserve linear features. With 
+<b>nearest</b> or <b>bilinear</b>, linear features can become zigzag 
+features after reprojection.
+
+<h2>KNOWN ISSUES</h2>
+Option <b>extent</b>=<em>region</em> works only when dataset has different
+projection than current location (<a href="r.proj.html">r.proj</a> is then used).
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.in.gdal.html">r.in.gdal</a>,<br>
+<a href="r.proj.html">r.proj</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Markus Metz<br>
+
+<p><i>Last changed: $Date: 2015-01-20 20:52:27 +0100 (Tue, 20 Jan 2015) $</i>

scripts/r.import/r.import.py

@@ -0,0 +1,335 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE:       r.import
+#
+# AUTHOR(S):    Markus Metz
+#
+# PURPOSE:      Import and reproject on the fly
+#
+# COPYRIGHT:    (C) 2015 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: Imports raster data into a GRASS raster map using GDAL library and reproject on the fly.
+#% keyword: raster
+#% keyword: import
+#% keyword: projection
+#%end
+#%option
+#% key: input
+#% type: string
+#% required: yes
+#% description: Name of GDAL dataset to be imported
+#% guisection: Input
+#%end
+#%option
+#% key: band
+#% type: integer
+#% required: no
+#% multiple: yes
+#% description: Input band(s) to select (default is all bands)
+#% guisection: Input
+#%end
+#%option
+#% key: memory
+#% type: integer
+#% required: no
+#% multiple: no
+#% options: 0-2047
+#% label: Maximum memory to be used (in MB)
+#% description: Cache size for raster rows
+#% answer: 300
+#%end
+#%option G_OPT_R_OUTPUT
+#% description: Name for output raster map
+#% guisection: Output
+#%end
+#%option
+#% key: resample
+#% type: string
+#% required: yes
+#% multiple: no
+#% options: nearest,bilinear,bicubic,lanczos,bilinear_f,bicubic_f,lanczos_f
+#% description: Resampling method to use for reprojection
+#% descriptions: nearest;nearest neighbor;bilinear;bilinear interpolation;bicubic;bicubic interpolation;lanczos;lanczos filter;bilinear_f;bilinear interpolation with fallback;bicubic_f;bicubic interpolation with fallback;lanczos_f;lanczos filter with fallback
+#% guisection: Output
+#%end
+#%option
+#% key: extent
+#% type: string
+#% required: no
+#% multiple: no
+#% options: input,region
+#% answer: input
+#% description: Output raster map extent
+#% descriptions: region;extent of current region;input;extent of input map
+#% guisection: Output
+#%end
+#%option
+#% key: resolution
+#% type: string
+#% required: no
+#% multiple: no
+#% answer: estimated
+#% options: estimated,value,region
+#% description: Resolution of output raster map (default: estimated)
+#% descriptions: estimated;estimated resolution;value;user-specified resolution;region;current region resolution
+#% guisection: Output
+#%end
+#%option
+#% key: resolution_value
+#% type: double
+#% required: no
+#% multiple: no
+#% description: Resolution of output raster map (use with option resolution=value)
+#% guisection: Output
+#%end
+#%flag
+#% key: e
+#% description: Estimate resolution only
+#% guisection: Optional
+#%end
+#%flag
+#% key: n
+#% description: Do not perform region cropping optimization
+#% guisection: Optional
+#%end
+
+
+import sys
+import os
+import atexit
+import math
+
+import grass.script as grass
+from grass.exceptions import CalledModuleError
+
+
+# initialize global vars
+TMPLOC = None
+SRCGISRC = None
+GISDBASE = None
+
+
+def cleanup():
+    # remove temp location
+    if TMPLOC:
+        grass.try_rmdir(os.path.join(GISDBASE, TMPLOC))
+    if SRCGISRC:
+        grass.try_remove(SRCGISRC)
+
+def main():
+    global TMPLOC, SRCGISRC, GISDBASE
+
+    GDALdatasource = options['input']
+    output = options['output']
+    method = options['resample']
+    memory = options['memory']
+    bands = options['band']
+    tgtres = options['resolution']
+    if options['resolution_value']:
+        if tgtres != 'value':
+            grass.fatal(_("To set custom resolution value, select 'value' in resolution option"))
+        tgtres_value = float(options['resolution_value'])
+        if tgtres_value <= 0:
+            grass.fatal(_("Resolution value can't be smaller than 0"))
+    elif tgtres == 'value':
+         grass.fatal(_("Please provide the resolution for the imported dataset or change to 'estimated' resolution"))
+
+    grassenv = grass.gisenv()
+    tgtloc = grassenv['LOCATION_NAME']
+    tgtmapset = grassenv['MAPSET']
+    GISDBASE = grassenv['GISDBASE']
+    tgtgisrc = os.environ['GISRC']
+    SRCGISRC = grass.tempfile()
+
+    TMPLOC = 'temp_import_location_' + str(os.getpid())
+
+    f = open(SRCGISRC, 'w')
+    f.write('MAPSET: PERMANENT\n')
+    f.write('GISDBASE: %s\n' % GISDBASE)
+    f.write('LOCATION_NAME: %s\n' % TMPLOC)
+    f.write('GUI: text\n')
+    f.close()
+
+    tgtsrs = grass.read_command('g.proj', flags='j', quiet=True)
+
+    # create temp location from input without import
+    grass.message(_("Creating temporary location for <%s>...") % GDALdatasource)
+    parameters = dict(input=GDALdatasource, output=output,
+                      memory=memory, flags='c',
+                      location=TMPLOC, quiet=True)
+    if bands:
+        parameters['band'] = bands
+    try:
+        grass.run_command('r.in.gdal', **parameters)
+    except CalledModuleError:
+        grass.fatal(_("Unable to read GDAL dataset <%s>") % GDALdatasource)
+
+    # switch to temp location
+    os.environ['GISRC'] = str(SRCGISRC)
+
+    # compare source and target srs
+    insrs = grass.read_command('g.proj', flags='j', quiet=True)
+
+    # switch to target location
+    os.environ['GISRC'] = str(tgtgisrc)
+
+    if insrs == tgtsrs:
+        # try r.in.gdal directly
+        parameters = dict(input=GDALdatasource, output=output,
+                          memory=memory, flags='k')
+        if bands:
+            parameters['band'] = bands
+        grass.message(_("Importing <%s>...") % GDALdatasource)
+        try:
+            grass.run_command('r.in.gdal', **parameters)
+        except CalledModuleError:
+            grass.fatal(_("Unable to import GDAL dataset <%s>") % GDALdatasource)
+        grass.message(_("Input <%s> successfully imported without reprojection") % GDALdatasource)
+        return 0
+
+    # make sure target is not xy
+    if grass.parse_command('g.proj', flags='g')['name'] == 'xy_location_unprojected':
+        grass.fatal(_("Coordinate reference system not available for current location <%s>") % tgtloc)
+
+    # switch to temp location
+    os.environ['GISRC'] = str(SRCGISRC)
+
+    # make sure input is not xy
+    if grass.parse_command('g.proj', flags='g')['name'] == 'xy_location_unprojected':
+        grass.fatal(_("Coordinate reference system not available for input <%s>") % GDALdatasource)
+
+    # import into temp location
+    grass.message(_("Importing <%s> to temporary location...") % GDALdatasource)
+    parameters = dict(input=GDALdatasource, output=output,
+                      memory=memory, flags='k')
+    if bands:
+        parameters['band'] = bands
+    try:
+        grass.run_command('r.in.gdal', **parameters)
+    except CalledModuleError:
+        grass.fatal(_("Unable to import GDAL dataset <%s>") % GDALdatasource)
+
+    outfiles = grass.list_grouped('raster')['PERMANENT']
+
+    # is output a group?
+    group = False
+    path = os.path.join(GISDBASE, TMPLOC, 'group', output)
+    if os.path.exists(path):
+        group = True
+        path = os.path.join(GISDBASE, TMPLOC, 'group', output, 'POINTS')
+        if os.path.exists(path):
+            grass.fatal(_("Input contains GCPs, rectification is required"))
+
+    # switch to target location
+    os.environ['GISRC'] = str(tgtgisrc)
+
+    region = grass.region()
+
+    rflags = None
+    if flags['n']:
+        rflags = 'n'
+
+    for outfile in outfiles:
+
+        n = region['n']
+        s = region['s']
+        e = region['e']
+        w = region['w']
+
+        grass.use_temp_region()
+
+        if options['extent'] == 'input':
+            # r.proj -g
+            try:
+                tgtextents = grass.read_command('r.proj', location=TMPLOC,
+                                                mapset='PERMANENT',
+                                                input=outfile, flags='g',
+                                                memory=memory, quiet=True)
+            except CalledModuleError:
+                grass.fatal(_("Unable to get reprojected map extent"))
+
+            srcregion = grass.parse_key_val(tgtextents, val_type=float, vsep=' ')
+            n = srcregion['n']
+            s = srcregion['s']
+            e = srcregion['e']
+            w = srcregion['w']
+
+            grass.run_command('g.region', n=n, s=s, e=e, w=w)
+
+        # v.in.region in tgt
+        vreg = 'vreg_' + str(os.getpid())
+        grass.run_command('v.in.region', output=vreg, quiet=True)
+
+        grass.del_temp_region()
+
+        # reproject to src
+        # switch to temp location
+        os.environ['GISRC'] = str(SRCGISRC)
+        try:
+            grass.run_command('v.proj', input=vreg, output=vreg,
+                              location=tgtloc, mapset=tgtmapset, quiet=True)
+        except CalledModuleError:
+            grass.fatal(_("Unable to reproject to source location"))
+
+        # set region from region vector
+        grass.run_command('g.region', raster=outfile)
+        grass.run_command('g.region', vector=vreg)
+        # align to fist band
+        grass.run_command('g.region', align=outfile)
+        # get number of cells
+        cells = grass.region()['cells']
+
+        estres = math.sqrt((n - s) * (e - w) / cells)
+
+        os.environ['GISRC'] = str(tgtgisrc)
+        grass.run_command('g.remove', type='vector', name=vreg,
+                          flags='f', quiet=True)
+
+        grass.message(_("Estimated target resolution for input band <%s>: %g") % (outfile, estres))
+        if flags['e']:
+            return 0
+
+        if options['extent'] == 'input':
+            grass.use_temp_region()
+            grass.run_command('g.region', n=n, s=s, e=e, w=w)
+
+        res = None
+        if tgtres == 'estimated':
+            res = estres
+        elif tgtres == 'value':
+            res = tgtres_value
+
+        # r.proj
+        grass.message(_("Reprojecting <%s>...") % outfile)
+        try:
+            grass.run_command('r.proj', location=TMPLOC,
+                              mapset='PERMANENT', input=outfile,
+                              method=method, resolution=res,
+                              memory=memory, flags=rflags, quiet=True)
+        except CalledModuleError:
+            grass.fatal(_("Unable to to reproject raster <%s>") % outfile)
+
+        if grass.raster_info(outfile)['min'] is None:
+            grass.fatal(_("The reprojected raster <%s> is empty") % outfile)
+
+        if options['extent'] == 'input':
+            grass.del_temp_region()
+
+    if group:
+        grass.run_command('i.group', group=output, input=','.join(outfiles))
+
+    return 0
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    atexit.register(cleanup)
+    sys.exit(main())