[GRASS-SVN] r62488 - grass-addons/grass7/raster/r.bioclim

svn_grass at osgeo.org svn_grass at osgeo.org
Thu Oct 30 03:12:34 PDT 2014


Author: mmetz
Date: 2014-10-30 03:12:34 -0700 (Thu, 30 Oct 2014)
New Revision: 62488

Modified:
   grass-addons/grass7/raster/r.bioclim/r.bioclim.html
   grass-addons/grass7/raster/r.bioclim/r.bioclim.py
Log:
r.bioclim: add option number of quartals

Modified: grass-addons/grass7/raster/r.bioclim/r.bioclim.html
===================================================================
--- grass-addons/grass7/raster/r.bioclim/r.bioclim.html	2014-10-30 10:10:46 UTC (rev 62487)
+++ grass-addons/grass7/raster/r.bioclim/r.bioclim.html	2014-10-30 10:12:34 UTC (rev 62488)
@@ -17,16 +17,17 @@
 <h4>Quarter years</h4>
 The bioclimatic indices referring to the wettest, driest, warmest or 
 coldest quarter differ from the official Worldclim bioclimatic indices 
-because <em><b>r.bioclim</b></em> divides a year into four quarters, 
-whereas Worldclim divides a year into 12 quarters. In 
-<em><b>r.bioclim</b></em>, the first quarter refers to Jan - Mar, the 
-second quarter to Apr - Jun, and the last quarter to Oct - Dec. In 
-Worldclim, the first quarter refers to Jan - Mar, the second quarter to 
-Feb - Apr, and the last quarter to Dec - Feb. The reason is that 
-Worldclim assumes long-term averages as input, whereas 
-<em><b>r.bioclim</b></em> assumes long-term-averages or records for a 
-specific year, in which case it does not make sense to calculate a 
-quarter by combining December, January, and February.
+if <em><b>r.bioclim</b></em> divides a year into four quarters 
+(<b>quarters=4</b>), whereas Worldclim divides a year into 12 quarters. 
+With <b>quarters=4</b>, the first quarter refers to Jan - Mar, the 
+second quarter to Apr - Jun, and the last quarter to Oct - Dec. 
+With <b>quarters=12</b>, the first quarter refers to Jan - Mar, 
+the second quarter to Feb - Apr, and the last quarter to Dec - Feb, as 
+in Worldclim. Using <b>quarters=12</b> makes thus only sense if 
+long-term averages are used as input. The default <b>quarters=4</b> 
+should be used with records for a specific year, because in this case 
+it does not make sense to calculate a quarter by combining December, 
+January, and February.
 
 <h4>List of bioclimatic indices</h4>
 

Modified: grass-addons/grass7/raster/r.bioclim/r.bioclim.py
===================================================================
--- grass-addons/grass7/raster/r.bioclim/r.bioclim.py	2014-10-30 10:10:46 UTC (rev 62487)
+++ grass-addons/grass7/raster/r.bioclim/r.bioclim.py	2014-10-30 10:12:34 UTC (rev 62488)
@@ -61,6 +61,13 @@
 #% description: Scale factor for output temperature
 #% answer: 10
 #%end
+#%option
+#% key: quartals
+#% type: integer
+#% description: Number of quartals to use 
+#% options: 4,12
+#% answer: 4
+#%end
 #%Option
 #% key: workers
 #% type: integer
@@ -78,8 +85,9 @@
 import atexit
 
 def cleanup():
-    if tmp:
-        grass.run_command('g.remove', flags = 'f', type = 'rast', pattern = tmp, quiet = True)
+    if tmp_pattern:
+        grass.run_command('g.remove', flags = 'f', type = 'rast', 
+	                  pattern = tmp_pattern, quiet = True)
 
 def main():
     tavg = options['tavg']
@@ -90,7 +98,10 @@
     tinscale = options['tinscale']
     toutscale = options['toutscale']
     workers = int(options['workers'])
+    quartals = int(options['quartals'])
 
+    qstep = 12 / quartals
+
     mapset = grass.gisenv()['MAPSET']
 
     # count input maps
@@ -114,7 +125,7 @@
     pid = os.getpid()
 
     # all temporary raster maps must follow this naming pattern
-    tmp = "%s.*.%d" % (outpre, pid)
+    tmp_pattern = "%s.*.%d" % (outpre, pid)
 
     tminl = tmin.split(',')
     tmaxl = tmax.split(',')
@@ -271,11 +282,18 @@
 
     # mean of mean for each quarter year
     grass.message(_("Mean temperature for each quarter year ..."))
-    for i in range(0, 12, 3):
-        tavgq = "%s.tavgq.%d.%d" % (outpre, i / 3, pid)
+    for i in range(quartals):
+        tavgq = "%s.tavgq.%02d.%d" % (outpre, i, pid)
 
+	m1 = i * qstep
+	m2 = m1 + 1
+	if m2 > 11:
+	    m2 = m2 - 12
+	m3 = m1 + 2
+	if m3 > 11:
+	    m3 = m3 - 12
         grass.run_command('r.series', 
-                           input = "%s,%s,%s" % (tavgl[i], tavgl[i + 1], tavgl[i + 2]),
+                           input = "%s,%s,%s" % (tavgl[m1], tavgl[m2], tavgl[m3]),
                            output = tavgq, method = 'average')
 
     # BIO10 = Mean Temperature of Warmest Quarter
@@ -286,7 +304,7 @@
     tavgq = grass.read_command("g.list",
                              quiet = True,
                              type = 'rast',
-                             pattern = '%s.tavgq.?.%d' % (outpre, pid),
+                             pattern = '%s.tavgq.??.%d' % (outpre, pid),
                              separator = ',',
                              mapset= '.')
 
@@ -302,40 +320,47 @@
                   oscale = toutscale,
                   biotmp = bio10,
                   iscale = tinscale)
+    grass.run_command('r.support', map = outpre + '.bio10',
+                      description = 'BIOCLIM10: Generated by r.bioclim')
+    grass.run_command('r.support', map = outpre + '.bio10',history = os.environ['CMDLINE'])
     grass.mapcalc("$bio = round(double($oscale) * $biotmp / $iscale)",
                   bio = outpre + '.bio11',
                   oscale = toutscale,
                   biotmp = bio11,
                   iscale = tinscale)
-    grass.run_command('r.support', map = outpre + '.bio10',
-                      description = 'BIOCLIM10: Generated by r.bioclim')
-    grass.run_command('r.support', map = outpre + '.bio10',history = os.environ['CMDLINE'])
     grass.run_command('r.support', map = outpre + '.bio11',
                       description = 'BIOCLIM11: Generated by r.bioclim')
     grass.run_command('r.support', map = outpre + '.bio11',history = os.environ['CMDLINE'])
+
     grass.run_command('g.remove', flags = 'f', type = 'rast', name = "%s,%s" % (bio10, bio11), quiet = True);
-    
 
     if not prec:
-        grass.run_command('g.remove', flags = 'f', type = 'rast', name = tmp, quiet = True)
+        grass.run_command('g.remove', flags = 'f', type = 'rast', 
+	                  pattern = tmp_pattern, quiet = True)
         sys.exit(1)
 
     precl = prec.split(',')
 
     # sum for each quarter year
     grass.message(_("Precipitation for each quarter year ..."))
-    for i in range(4):
-        precq = "%s.precq.%d.%d" % (outpre, i + 1, pid)
+    for i in range(quartals):
+        precq = "%s.precq.%02d.%d" % (outpre, i + 1, pid)
 
-        j = i * 3
+        m1 = i * qstep
+	m2 = m1 + 1
+	if m2 > 11:
+	    m2 = m2 - 12
+	m3 = m1 + 2
+	if m3 > 11:
+	    m3 = m3 - 12
         grass.run_command('r.series', 
-                           input = "%s,%s,%s" % (precl[j], precl[j + 1], precl[j + 2]),
+                           input = "%s,%s,%s" % (precl[m1], precl[m2], precl[m3]),
                            output = precq, method = 'sum')
 
     precq = grass.read_command("g.list",
                              quiet = True,
                              type = 'rast',
-                             pattern = '%s.precq.?.%d' % (outpre, pid),
+                             pattern = '%s.precq.??.%d' % (outpre, pid),
                              separator = ',',
                              mapset= '.')
 
@@ -361,16 +386,42 @@
     # BIO8 = Mean Temperature of Wettest Quarter
     grass.message(_("BIO8 = Mean Temperature of Wettest Quarter ..."))
 
-    grass.mapcalc("$bio = round(if($wettestq == 0, $tavgq0, \
-                                if($wettestq == 1, $tavgq1, \
-                                if($wettestq == 2, $tavgq2, \
-                                if($wettestq == 3, $tavgq3, null())))) \
-				* $oscale / $iscale)",
-                  bio = outpre + '.bio08',
-                  wettestq = wettestq, 
-                  tavgq0 = tavgql[0], tavgq1 = tavgql[1],
-                  tavgq2 = tavgql[2], tavgq3 = tavgql[3],
-                  oscale = toutscale, iscale = tinscale)
+    if quartals == 4:
+	grass.mapcalc("$bio = round(if($wettestq == 0, $tavgq0, \
+				    if($wettestq == 1, $tavgq1, \
+				    if($wettestq == 2, $tavgq2, \
+				    if($wettestq == 3, $tavgq3, null())))) \
+				    * $oscale / $iscale)",
+		      bio = outpre + '.bio08',
+		      wettestq = wettestq, 
+		      tavgq0 = tavgql[0], tavgq1 = tavgql[1],
+		      tavgq2 = tavgql[2], tavgq3 = tavgql[3],
+		      oscale = toutscale, iscale = tinscale)
+    else:
+	# quartals == 12
+	grass.mapcalc("$bio = round(if($wettestq ==  0, $tavgq0, \
+				    if($wettestq ==  1, $tavgq1, \
+				    if($wettestq ==  2, $tavgq2, \
+				    if($wettestq ==  3, $tavgq3, \
+				    if($wettestq ==  4, $tavgq4, \
+				    if($wettestq ==  5, $tavgq5, \
+				    if($wettestq ==  6, $tavgq6, \
+				    if($wettestq ==  7, $tavgq7, \
+				    if($wettestq ==  8, $tavgq8, \
+				    if($wettestq ==  9, $tavgq9, \
+				    if($wettestq == 10, $tavgq10, \
+				    if($wettestq == 11, $tavgq11, null())))))))))))) \
+				    * $oscale / $iscale)",
+		      bio = outpre + '.bio08',
+		      wettestq = wettestq, 
+		      tavgq0 = tavgql[0], tavgq1 = tavgql[1],
+		      tavgq2 = tavgql[2], tavgq3 = tavgql[3],
+		      tavgq4 = tavgql[4], tavgq5 = tavgql[5],
+		      tavgq6 = tavgql[6], tavgq7 = tavgql[7],
+		      tavgq8 = tavgql[8], tavgq9 = tavgql[9],
+		      tavgq10 = tavgql[10], tavgq11 = tavgql[11],
+		      oscale = toutscale, iscale = tinscale)
+
     grass.run_command('r.support', map = outpre + '.bio08',
                       description = 'BIOCLIM08: Generated by r.bioclim')
     grass.run_command('r.support', map = outpre + '.bio08',history = os.environ['CMDLINE'])
@@ -378,16 +429,42 @@
     # BIO9 = Mean Temperature of Driest Quarter
     grass.message(_("BIO9 = Mean Temperature of Driest Quarter ..."))
 
-    grass.mapcalc("$bio = round(if($driestq == 0, $tavgq0, \
-                                if($driestq == 1, $tavgq1, \
-                                if($driestq == 2, $tavgq2, \
-                                if($driestq == 3, $tavgq3, null())))) \
-				* $oscale / $iscale)",
-                  bio = outpre + '.bio09',
-                  driestq = driestq, 
-                  tavgq0 = tavgql[0], tavgq1 = tavgql[1],
-                  tavgq2 = tavgql[2], tavgq3 = tavgql[3],
-                  oscale = toutscale, iscale = tinscale)
+    if quartals == 4:
+	grass.mapcalc("$bio = round(if($driestq == 0, $tavgq0, \
+				    if($driestq == 1, $tavgq1, \
+				    if($driestq == 2, $tavgq2, \
+				    if($driestq == 3, $tavgq3, null())))) \
+				    * $oscale / $iscale)",
+		      bio = outpre + '.bio09',
+		      driestq = driestq, 
+		      tavgq0 = tavgql[0], tavgq1 = tavgql[1],
+		      tavgq2 = tavgql[2], tavgq3 = tavgql[3],
+		      oscale = toutscale, iscale = tinscale)
+    else:
+	# quartals == 12
+	grass.mapcalc("$bio = round(if($driestq ==  0, $tavgq0, \
+				    if($driestq ==  1, $tavgq1, \
+				    if($driestq ==  2, $tavgq2, \
+				    if($driestq ==  3, $tavgq3, \
+				    if($driestq ==  4, $tavgq4, \
+				    if($driestq ==  5, $tavgq5, \
+				    if($driestq ==  6, $tavgq6, \
+				    if($driestq ==  7, $tavgq7, \
+				    if($driestq ==  8, $tavgq8, \
+				    if($driestq ==  9, $tavgq9, \
+				    if($driestq == 10, $tavgq10, \
+				    if($driestq == 11, $tavgq11, null())))))))))))) \
+				    * $oscale / $iscale)",
+		      bio = outpre + '.bio09',
+		      driestq = driestq, 
+		      tavgq0 = tavgql[0], tavgq1 = tavgql[1],
+		      tavgq2 = tavgql[2], tavgq3 = tavgql[3],
+		      tavgq4 = tavgql[4], tavgq5 = tavgql[5],
+		      tavgq6 = tavgql[6], tavgq7 = tavgql[7],
+		      tavgq8 = tavgql[8], tavgq9 = tavgql[9],
+		      tavgq10 = tavgql[10], tavgq11 = tavgql[11],
+		      oscale = toutscale, iscale = tinscale)
+
     grass.run_command('r.support', map = outpre + '.bio09',
                       description = 'BIOCLIM09: Generated by r.bioclim')
     grass.run_command('r.support', map = outpre + '.bio09',history = os.environ['CMDLINE'])
@@ -455,14 +532,38 @@
     # BIO18 = Precipitation of Warmest Quarter
     grass.message(_("BIO18 = Precipitation of Warmest Quarter ..."))
 
-    grass.mapcalc("$bio = round(if($warmestq == 0, $precq0, \
-                                if($warmestq == 1, $precq1, \
-                                if($warmestq == 2, $precq2, \
-                                if($warmestq == 3, $precq3, null())))))",
-                  bio = outpre + '.bio18',
-                  warmestq = warmestq, 
-                  precq0 = precql[0], precq1 = precql[1],
-                  precq2 = precql[2], precq3 = precql[3])
+    if quartals == 4:
+	grass.mapcalc("$bio = round(if($warmestq == 0, $precq0, \
+				    if($warmestq == 1, $precq1, \
+				    if($warmestq == 2, $precq2, \
+				    if($warmestq == 3, $precq3, null())))))",
+		      bio = outpre + '.bio18',
+		      warmestq = warmestq, 
+		      precq0 = precql[0], precq1 = precql[1],
+		      precq2 = precql[2], precq3 = precql[3])
+    else:
+	# quartals == 12
+	grass.mapcalc("$bio = round(if($warmestq ==  0, $precq0, \
+				    if($warmestq ==  1, $precq1, \
+				    if($warmestq ==  2, $precq2, \
+				    if($warmestq ==  3, $precq3, \
+				    if($warmestq ==  4, $precq4, \
+				    if($warmestq ==  5, $precq5, \
+				    if($warmestq ==  6, $precq6, \
+				    if($warmestq ==  7, $precq7, \
+				    if($warmestq ==  8, $precq8, \
+				    if($warmestq ==  9, $precq9, \
+				    if($warmestq == 10, $precq10, \
+				    if($warmestq == 11, $precq11, null())))))))))))))",
+		      bio = outpre + '.bio18',
+		      warmestq = warmestq, 
+		      precq0 = precql[0], precq1 = precql[1],
+		      precq2 = precql[2], precq3 = precql[3],
+		      precq4 = precql[4], precq5 = precql[5],
+		      precq6 = precql[6], precq7 = precql[7],
+		      precq8 = precql[8], precq9 = precql[9],
+		      precq10 = precql[10], precq11 = precql[11])
+
     grass.run_command('r.support', map = outpre + '.bio18',
                       description = 'BIOCLIM18: Generated by r.bioclim')
     grass.run_command('r.support', map = outpre + '.bio18',history = os.environ['CMDLINE'])
@@ -470,23 +571,48 @@
     # BIO19 = Precipitation of Coldest Quarter
     grass.message(_("BIO19 = Precipitation of Coldest Quarter ..."))
 
-    grass.mapcalc("$bio = round(if($coldestq == 0, $precq0, \
-                                if($coldestq == 1, $precq1, \
-                                if($coldestq == 2, $precq2, \
-                                if($coldestq == 3, $precq3, null())))))",
-                  bio = outpre + '.bio19',
-                  coldestq = coldestq, 
-                  precq0 = precql[0], precq1 = precql[1],
-                  precq2 = precql[2], precq3 = precql[3])
+    if quartals == 4:
+	grass.mapcalc("$bio = round(if($coldestq == 0, $precq0, \
+				    if($coldestq == 1, $precq1, \
+				    if($coldestq == 2, $precq2, \
+				    if($coldestq == 3, $precq3, null())))))",
+		      bio = outpre + '.bio19',
+		      coldestq = coldestq, 
+		      precq0 = precql[0], precq1 = precql[1],
+		      precq2 = precql[2], precq3 = precql[3])
+    else:
+	# quartals == 12
+	grass.mapcalc("$bio = round(if($coldestq ==  0, $precq0, \
+				    if($coldestq ==  1, $precq1, \
+				    if($coldestq ==  2, $precq2, \
+				    if($coldestq ==  3, $precq3, \
+				    if($coldestq ==  4, $precq4, \
+				    if($coldestq ==  5, $precq5, \
+				    if($coldestq ==  6, $precq6, \
+				    if($coldestq ==  7, $precq7, \
+				    if($coldestq ==  8, $precq8, \
+				    if($coldestq ==  9, $precq9, \
+				    if($coldestq == 10, $precq10, \
+				    if($coldestq == 11, $precq11, null())))))))))))))",
+		      bio = outpre + '.bio19',
+		      coldestq = coldestq, 
+		      precq0 = precql[0], precq1 = precql[1],
+		      precq2 = precql[2], precq3 = precql[3],
+		      precq4 = precql[4], precq5 = precql[5],
+		      precq6 = precql[6], precq7 = precql[7],
+		      precq8 = precql[8], precq9 = precql[9],
+		      precq10 = precql[10], precq11 = precql[11])
+
     grass.run_command('r.support', map = outpre + '.bio19',
                       description = 'BIOCLIM19: Generated by r.bioclim')
     grass.run_command('r.support', map = outpre + '.bio19',history = os.environ['CMDLINE'])
 
-    grass.run_command('g.remove', flags = 'f', type = 'rast', name = tmp, quiet = True)
+    grass.run_command('g.remove', flags = 'f', type = 'rast', 
+                      pattern = tmp_pattern, quiet = True)
 
 
 if __name__ == "__main__":
     options, flags = grass.parser()
-    tmp = None
+    tmp_pattern = None
     atexit.register(cleanup)
     main()



More information about the grass-commit mailing list