[GRASS-SVN] r33649 - grass/trunk/scripts/i.in.spotvgt
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Oct 1 20:32:11 EDT 2008
Author: glynn
Date: 2008-10-01 20:32:10 -0400 (Wed, 01 Oct 2008)
New Revision: 33649
Added:
grass/trunk/scripts/i.in.spotvgt/i.in.spotvgt.py
Log:
Convert i.in.spotvgt to Python (untested)
Added: grass/trunk/scripts/i.in.spotvgt/i.in.spotvgt.py
===================================================================
--- grass/trunk/scripts/i.in.spotvgt/i.in.spotvgt.py (rev 0)
+++ grass/trunk/scripts/i.in.spotvgt/i.in.spotvgt.py 2008-10-02 00:32:10 UTC (rev 33649)
@@ -0,0 +1,271 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: i.in.spot
+#
+# AUTHOR(S): Markus Neteler
+# Converted to Python by Glynn Clements
+#
+# PURPOSE: Import SPOT VEGETATION data into a GRASS raster map
+# SPOT Vegetation (1km, global) data:
+# http://free.vgt.vito.be/
+#
+# COPYRIGHT: (c) 2004-2008 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.
+#
+#############################################################################
+#
+# REQUIREMENTS:
+# - gdal: http://www.gdal.org
+#
+# Notes:
+# * According to the faq (http://www.vgt.vito.be/faq/faq.html), SPOT vegetation
+# coordinates refer to the center of a pixel.
+# * GDAL coordinates refer to the corner of a pixel
+# -> correction of 0001/0001_LOG.TXT coordinates by 0.5 pixel
+#############################################################################
+
+#%Module
+#% description: Import of SPOT VGT NDVI file into a raster map
+#% keywords: raster, imagery, import
+#%End
+#%flag
+#% key: a
+#% description: also import quality map (SM status map layer) and filter NDVI map
+#%end
+#%option
+#% key: file
+#% type: string
+#% description: existing SPOT VGT NDVI HDF file (0001_NDV.HDF)
+#% gisprompt: old_file,file,input
+#% required : yes
+#%end
+#%option
+#% key: rast
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Name for output raster map
+#% required : no
+#%end
+
+import sys
+import os
+import subprocess
+import atexit
+import grass
+
+vrt = """<VRTDataset rasterXSize="$XSIZE" rasterYSize="$YSIZE">
+ <SRS>GEOGCS["wgs84",DATUM["WGS_1984",SPHEROID["wgs84",6378137,298.257223563],TOWGS84[0.000,0.000,0.000]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]</SRS>
+ <GeoTransform>$WESTCORNER, $RESOLUTION, 0.0000000000000000e+00, $NORTHCORNER, 0.0000000000000e+00, -$RESOLUTION</GeoTransform>
+ <VRTRasterBand dataType="Byte" band="1">
+ <NoDataValue>0.00000000000000E+00</NoDataValue>
+ <ColorInterp>Gray</ColorInterp>
+ <SimpleSource>
+ <SourceFilename>$FILENAME</SourceFilename>
+ <SourceBand>1</SourceBand>
+ <SrcRect xOff="0" yOff="0" xSize="$XSIZE" ySize="$YSIZE"/>
+ <DstRect xOff="0" yOff="0" xSize="$XSIZE" ySize="$YSIZE"/>
+ </SimpleSource>
+ </VRTRasterBand>
+</VRTDataset>"""
+
+#### a function for writing VRT files
+def create_VRT_file(projfile, vrtfile, infile):
+ fh = file(projfile)
+ kv = {}
+ for l in fh:
+ f = l.rstrip('\r\n').split()
+ if f < 2:
+ continue
+ kv[f[0]] = f[1]
+ fh.close()
+
+ north_center = kv['CARTO_UPPER_LEFT_Y']
+ south_center = kv['CARTO_LOWER_LEFT_Y']
+ east_center = kv['CARTO_UPPER_RIGHT_X']
+ west_center = kv['CARTO_UPPER_LEFT_X']
+ map_proj_res = kv['MAP_PROJ_RESOLUTION']
+ xsize = kv['IMAGE_UPPER_RIGHT_COL']
+ ysize = kv['IMAGE_LOWER_RIGHT_ROW']
+
+ resolution = float(map_proj_res)
+ north_corner = float(north_center) + resolution / 2
+ south_corner = float(south_center) - resolution / 2
+ east_corner = float(east_center ) + resolution / 2
+ west_corner = float(west_center ) - resolution / 2
+
+ t = string.Template(vrt)
+ s = t.substitute(NORTHCORNER = north_corner, WESTCORNER = west_corner,
+ XSIZE = xsize, YSIZE = ysize, RESOLUTION = map_proj_res,
+ FILENAME = infile)
+ outf = file(vrtfile, 'w')
+ outf.write(s)
+ outf.close()
+
+def cleanup():
+ #### clean up the mess
+ grass.try_remove(vrtfile)
+ grass.try_remove(tmpfile)
+
+def main():
+ global vrtfile, tmpfile
+
+ infile = options['file']
+ rast = options['rast']
+ also = flags['a']
+
+ #### check for gdalinfo (just to check if installation is complete)
+ if not grass.find_program('gdalinfo', ['--version']):
+ grass.fatal("'gdalinfo' not found, install GDAL tools first (http://www.gdal.org)")
+
+ pid = str(os.getpid())
+ tmpfile = grass.tempfile()
+
+ ################### let's go
+
+ spotdir = os.path.dirname(infile)
+ spotname = grass.basename(infile, 'hdf')
+
+ if rast:
+ name = rast
+ else:
+ name = spotname
+
+ if not grass.overwrite() and grass.find_file(name)['file']:
+ grass.fatal("<%s> already exists. Aborting." % name)
+
+ # still a ZIP file? (is this portable?? see the r.in.srtm script for ideas)
+ if infile.lower().endswith('.zip'):
+ grass.fatal("Please extract %s before import." % infile)
+
+ try:
+ p = subprocess.Popen(['file', '-ib', infile], stdout = subprocess.PIPE)
+ s = p.communicate()[0]
+ if s == "application/x-zip":
+ grass.fatal("Please extract %s before import." % infile)
+ except:
+ pass
+
+ ### create VRT header for NDVI
+
+ projfile = os.path.join(spotdir, "0001_LOG.TXT")
+ vrtfile = tmpfile + '.vrt'
+
+ # first process the NDVI:
+ grass.try_remove(vrtfile)
+ create_VRT_file(projfile, vrtfile, infile)
+
+ ## let's import the NDVI map...
+ grass.message("Importing SPOT VGT NDVI map...")
+ if grass.run_command('r.in.gdal', input = vrtfile, output = name) != 0:
+ grass.fatal("An error occurred. Stop.")
+
+ grass.message("Imported SPOT VEGETATION NDVI map <%s>." % name)
+
+ #################
+ ## http://www.vgt.vito.be/faq/FAQS/faq19.html
+ # What is the relation between the digital number and the real NDVI ?
+ # Real NDVI =coefficient a * Digital Number + coefficient b
+ # = a * DN +b
+ #
+ # Coefficient a = 0.004
+ # Coefficient b = -0.1
+
+ # clone current region
+ # switch to a temporary region
+ grass.use_temp_region()
+
+ grass.run_command('g.region', rast = name, quiet = True)
+
+ grass.message("Remapping digital numbers to NDVI...")
+ tmpname = "%s_%s" % (name, pid)
+ e = "%s = 0.004 * %s - 0.1" % (tmpname, name)
+ grass.run_command('r.mapcalc', expression = e)
+ grass.run_command('g.remove', rast = name, quiet = True)
+ grass.run_command('g.rename', rast = (tmpname, name), quiet = True)
+
+ #apply color table:
+ grass.run_command('r.colors', map = name, color = 'ndvi', quiet = True)
+
+ ##########################
+ # second, optionally process the SM quality map:
+
+ #SM Status Map
+ # http://www.vgt.vito.be/faq/FAQS/faq15.html
+ #Data about
+ # Bit NR 7: Radiometric quality for B0 coded as 0 if bad and 1 if good
+ # Bit NR 6: Radiometric quality for B2 coded as 0 if bad and 1 if good
+ # Bit NR 5: Radiometric quality for B3 coded as 0 if bad and 1 if good
+ # Bit NR 4: Radiometric quality for MIR coded as 0 if bad and 1 if good
+ # Bit NR 3: land code 1 or water code 0
+ # Bit NR 2: ice/snow code 1 , code 0 if there is no ice/snow
+ # Bit NR 1: 0 0 1 1
+ # Bit NR 0: 0 1 0 1
+ # clear shadow uncertain cloud
+ #
+ #Note:
+ # pos 7 6 5 4 3 2 1 0 (bit position)
+ # 128 64 32 16 8 4 2 1 (values for 8 bit)
+ #
+ #
+ # Bit 4-7 should be 1: their sum is 240
+ # Bit 3 land code, should be 1, sum up to 248 along with higher bits
+ # Bit 2 ice/snow code
+ # Bit 0-1 should be 0
+ #
+ # A good map threshold: >= 248
+
+ if also:
+ grass.message("Importing SPOT VGT NDVI quality map...")
+ grass.try_remove(vrtfile)
+ qfile = infile.replace('NDV','SM')
+ create_VRT_file(projfile, vrtfile, qfile)
+
+ ## let's import the SM quality map...
+ smfile = name + '.sm'
+ if grass.run_command('r.in.gdal', input = vrtfile, output = smfile) != 0:
+ grass.fatal("An error occurred. Stop.")
+
+ # some of the possible values:
+ rules = [r + '\n' for r in [
+ '8 50 50 50',
+ '11 70 70 70',
+ '12 90 90 90',
+ '60 grey',
+ '155 blue',
+ '232 violet',
+ '235 red',
+ '236 brown',
+ '248 orange',
+ '251 yellow',
+ '252 green'
+ ]]
+ grass.write_command('r.colors', map = smfile, rules = '-', stdin = rules)
+
+ grass.message("Imported SPOT VEGETATION SM quality map <%s>." % smfile)
+ grass.message("Note: A snow map can be extracted by category 252 (d.rast %s cat=252)" % smfile)
+ grass.message("")
+ grass.message("Filtering NDVI map by Status Map quality layer...")
+
+ filtfile = "%s_filt" % name
+ e = "%s = if(%s >= 248, %s, null())" % (filtfile, smfile, name)
+ grass.run_command('r.mapcalc', expression = e)
+ grass.run_command('r.colors', map = filtfile, color = 'ndvi', quiet = True)
+ grass.message("Filtered SPOT VEGETATION NDVI map <%s>." % filtfile)
+
+ # write cmd history:
+ grass.raster_history(name)
+ grass.raster_history(smfile)
+ grass.raster_history(filtfile)
+
+ grass.message("Done.")
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ atexit.register(cleanup)
+ main()
+
Property changes on: grass/trunk/scripts/i.in.spotvgt/i.in.spotvgt.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list