1#!/usr/local/bin/python3.8 2 3############################################################################ 4# 5# MODULE: i.in.spot 6# 7# AUTHOR(S): Markus Neteler 8# Converted to Python by Glynn Clements 9# 10# PURPOSE: Import SPOT VEGETATION data into a GRASS raster map 11# SPOT Vegetation (1km, global) data: 12# http://free.vgt.vito.be/ 13# 14# COPYRIGHT: (c) 2004-2011 GRASS Development Team 15# 16# This program is free software under the GNU General Public 17# License (>=v2). Read the file COPYING that comes with GRASS 18# for details. 19# 20############################################################################# 21# 22# REQUIREMENTS: 23# - gdal: http://www.gdal.org 24# 25# Notes: 26# * According to the faq (http://www.vgt.vito.be/faq/faq.html), SPOT vegetation 27# coordinates refer to the center of a pixel. 28# * GDAL coordinates refer to the corner of a pixel 29# -> correction of 0001/0001_LOG.TXT coordinates by 0.5 pixel 30############################################################################# 31 32#%module 33#% description: Imports SPOT VGT NDVI data into a raster map. 34#% keyword: imagery 35#% keyword: import 36#% keyword: NDVI 37#% keyword: SPOT 38#%end 39#%flag 40#% key: a 41#% description: Also import quality map (SM status map layer) and filter NDVI map 42#%end 43#%option G_OPT_F_INPUT 44#% description: Name of input SPOT VGT NDVI HDF file 45#%end 46#% option G_OPT_R_OUTPUT 47#% required : no 48#%end 49 50import os 51import atexit 52import string 53import grass.script as gscript 54from grass.exceptions import CalledModuleError 55 56 57vrt = """<VRTDataset rasterXSize="$XSIZE" rasterYSize="$YSIZE"> 58 <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> 59 <GeoTransform>$WESTCORNER, $RESOLUTION, 0.0000000000000000e+00, $NORTHCORNER, 0.0000000000000e+00, -$RESOLUTION</GeoTransform> 60 <VRTRasterBand dataType="Byte" band="1"> 61 <NoDataValue>0.00000000000000E+00</NoDataValue> 62 <ColorInterp>Gray</ColorInterp> 63 <SimpleSource> 64 <SourceFilename>$FILENAME</SourceFilename> 65 <SourceBand>1</SourceBand> 66 <SrcRect xOff="0" yOff="0" xSize="$XSIZE" ySize="$YSIZE"/> 67 <DstRect xOff="0" yOff="0" xSize="$XSIZE" ySize="$YSIZE"/> 68 </SimpleSource> 69 </VRTRasterBand> 70</VRTDataset>""" 71 72# a function for writing VRT files 73 74 75def create_VRT_file(projfile, vrtfile, infile): 76 fh = open(projfile) 77 kv = {} 78 for l in fh: 79 f = l.rstrip('\r\n').split() 80 if f < 2: 81 continue 82 kv[f[0]] = f[1] 83 fh.close() 84 85 north_center = kv['CARTO_UPPER_LEFT_Y'] 86 # south_center = kv['CARTO_LOWER_LEFT_Y'] 87 # east_center = kv['CARTO_UPPER_RIGHT_X'] 88 west_center = kv['CARTO_UPPER_LEFT_X'] 89 map_proj_res = kv['MAP_PROJ_RESOLUTION'] 90 xsize = kv['IMAGE_UPPER_RIGHT_COL'] 91 ysize = kv['IMAGE_LOWER_RIGHT_ROW'] 92 93 resolution = float(map_proj_res) 94 north_corner = float(north_center) + resolution / 2 95 # south_corner = float(south_center) - resolution / 2 96 # east_corner = float(east_center) + resolution / 2 97 west_corner = float(west_center) - resolution / 2 98 99 t = string.Template(vrt) 100 s = t.substitute(NORTHCORNER=north_corner, WESTCORNER=west_corner, 101 XSIZE=xsize, YSIZE=ysize, RESOLUTION=map_proj_res, 102 FILENAME=infile) 103 outf = open(vrtfile, 'w') 104 outf.write(s) 105 outf.close() 106 107 108def cleanup(): 109 # clean up the mess 110 gscript.try_remove(vrtfile) 111 gscript.try_remove(tmpfile) 112 113 114def main(): 115 global vrtfile, tmpfile 116 117 infile = options['input'] 118 rast = options['output'] 119 also = flags['a'] 120 121 # check for gdalinfo (just to check if installation is complete) 122 if not gscript.find_program('gdalinfo', '--help'): 123 gscript.fatal(_("'gdalinfo' not found, install GDAL tools first " 124 "(http://www.gdal.org)")) 125 126 pid = str(os.getpid()) 127 tmpfile = gscript.tempfile() 128 129 # let's go 130 131 spotdir = os.path.dirname(infile) 132 spotname = gscript.basename(infile, 'hdf') 133 134 if rast: 135 name = rast 136 else: 137 name = spotname 138 139 if not gscript.overwrite() and gscript.find_file(name)['file']: 140 gscript.fatal(_("<%s> already exists. Aborting.") % name) 141 142 # still a ZIP file? (is this portable?? see the r.in.srtm script for 143 # ideas) 144 if infile.lower().endswith('.zip'): 145 gscript.fatal(_("Please extract %s before import.") % infile) 146 147 try: 148 p = gscript.Popen(['file', '-ib', infile], stdout=gscript.PIPE) 149 s = p.communicate()[0] 150 if s == "application/x-zip": 151 gscript.fatal(_("Please extract %s before import.") % infile) 152 except: 153 pass 154 155 # create VRT header for NDVI 156 157 projfile = os.path.join(spotdir, "0001_LOG.TXT") 158 vrtfile = tmpfile + '.vrt' 159 160 # first process the NDVI: 161 gscript.try_remove(vrtfile) 162 create_VRT_file(projfile, vrtfile, infile) 163 164 # let's import the NDVI map... 165 gscript.message(_("Importing SPOT VGT NDVI map...")) 166 try: 167 gscript.run_command('r.in.gdal', input=vrtfile, output=name) 168 except CalledModuleError: 169 gscript.fatal(_("An error occurred. Stop.")) 170 171 gscript.message(_("Imported SPOT VEGETATION NDVI map <%s>.") % name) 172 173 ################# 174 # http://www.vgt.vito.be/faq/FAQS/faq19.html 175 # What is the relation between the digital number and the real NDVI ? 176 # Real NDVI =coefficient a * Digital Number + coefficient b 177 # = a * DN +b 178 # 179 # Coefficient a = 0.004 180 # Coefficient b = -0.1 181 182 # clone current region 183 # switch to a temporary region 184 gscript.use_temp_region() 185 186 gscript.run_command('g.region', raster=name, quiet=True) 187 188 gscript.message(_("Remapping digital numbers to NDVI...")) 189 tmpname = "%s_%s" % (name, pid) 190 gscript.mapcalc("$tmpname = 0.004 * $name - 0.1", tmpname=tmpname, name=name) 191 gscript.run_command('g.remove', type='raster', name=name, quiet=True, 192 flags='f') 193 gscript.run_command('g.rename', raster=(tmpname, name), quiet=True) 194 195 # write cmd history: 196 gscript.raster_history(name) 197 198 # apply color table: 199 gscript.run_command('r.colors', map=name, color='ndvi', quiet=True) 200 201 ########################## 202 # second, optionally process the SM quality map: 203 204 # SM Status Map 205 # http://nieuw.vgt.vito.be/faq/FAQS/faq22.html 206 # Data about 207 # Bit NR 7: Radiometric quality for B0 coded as 0 if bad and 1 if good 208 # Bit NR 6: Radiometric quality for B2 coded as 0 if bad and 1 if good 209 # Bit NR 5: Radiometric quality for B3 coded as 0 if bad and 1 if good 210 # Bit NR 4: Radiometric quality for MIR coded as 0 if bad and 1 if good 211 # Bit NR 3: land code 1 or water code 0 212 # Bit NR 2: ice/snow code 1 , code 0 if there is no ice/snow 213 # Bit NR 1: 0 0 1 1 214 # Bit NR 0: 0 1 0 1 215 # clear shadow uncertain cloud 216 # 217 # Note: 218 # pos 7 6 5 4 3 2 1 0 (bit position) 219 # 128 64 32 16 8 4 2 1 (values for 8 bit) 220 # 221 # 222 # Bit 4-7 should be 1: their sum is 240 223 # Bit 3 land code, should be 1, sum up to 248 along with higher bits 224 # Bit 2 ice/snow code 225 # Bit 0-1 should be 0 226 # 227 # A good map threshold: >= 248 228 229 if also: 230 gscript.message(_("Importing SPOT VGT NDVI quality map...")) 231 gscript.try_remove(vrtfile) 232 qname = spotname.replace('NDV', 'SM') 233 qfile = os.path.join(spotdir, qname) 234 create_VRT_file(projfile, vrtfile, qfile) 235 236 # let's import the SM quality map... 237 smfile = name + '.sm' 238 try: 239 gscript.run_command('r.in.gdal', input=vrtfile, output=smfile) 240 except CalledModuleError: 241 gscript.fatal(_("An error occurred. Stop.")) 242 243 # some of the possible values: 244 rules = [r + '\n' for r in ['8 50 50 50', 245 '11 70 70 70', 246 '12 90 90 90', 247 '60 grey', 248 '155 blue', 249 '232 violet', 250 '235 red', 251 '236 brown', 252 '248 orange', 253 '251 yellow', 254 '252 green']] 255 gscript.write_command('r.colors', map=smfile, rules='-', stdin=rules) 256 257 gscript.message(_("Imported SPOT VEGETATION SM quality map <%s>.") % 258 smfile) 259 gscript.message(_("Note: A snow map can be extracted by category " 260 "252 (d.rast %s cat=252)") % smfile) 261 gscript.message("") 262 gscript.message(_("Filtering NDVI map by Status Map quality layer...")) 263 264 filtfile = "%s_filt" % name 265 gscript.mapcalc("$filtfile = if($smfile % 4 == 3 || " 266 "($smfile / 16) % 16 == 0, null(), $name)", 267 filtfile=filtfile, smfile=smfile, name=name) 268 gscript.run_command('r.colors', map=filtfile, color='ndvi', quiet=True) 269 gscript.message(_("Filtered SPOT VEGETATION NDVI map <%s>.") % 270 filtfile) 271 272 # write cmd history: 273 gscript.raster_history(smfile) 274 gscript.raster_history(filtfile) 275 276 gscript.message(_("Done.")) 277 278if __name__ == "__main__": 279 options, flags = gscript.parser() 280 atexit.register(cleanup) 281 main() 282