[GRASS-user] region issues using python multiprocessing.Pool()

[Eric Goddard]

Hi all,

I'm attempting to convert WorldView-2 DNs to top of atmosphere
reflectance using the formulas provided in
http://www.digitalglobe.com/sites/default/files/Radiometric_Use_of_WorldView-2_Imagery%20%281%29.pdf.

If I run the script without multiprocessing enabled the script
executes properly; If I run it with multiprocessing enabled, The
script executes successfully but the  pixel range is 0-0; I think this
is caused by the region settings. Can the use_temp_region() function
be used as it is below to execute multiple mapcalc jobs in parallel
that have different extents? Once I have this working I'd like to make
a more generic tool a-la the i.landsat/toar tools.

Thanks,
Eric

r.info from one of the parrallel rasters:
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    DCELL                                                      |
 |   Rows:         4096                                                       |
 |   Columns:      4096                                                       |
 |   Total Cells:  16777216                                                   |
 |        Projection: Lambert Conformal Conic                                 |
 |            N: 335872.03166667    S: 308995.445   Res: 6.56166667           |
 |            E:  798935.41    W: 772058.82333333   Res: 6.56166667           |
 |   Range of data:    min = 0  max = 0                                       |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.mapcalc                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    0.009295654 * x2NOV04172233_M3DS_R4C1_052823926030_01_P001.1 *          |
 |    0.99154845 ^ 2 * 3.1415927 / (1758.2229 * cos(50.9))

The functions:

def CalculateTOAR(rastItems):
   raster = rastItems[0]
   absCalFactor = rastItems[1]
   theta = rastItems[2]
   earthSunDist = rastItems[3]
   eSun = rastItems[4]
   output = "%s_TOAR%s" % (raster[:-2], raster[-2:])

   grass.use_temp_region()
   grass.run_command('g.region', rast=raster)

   grass.mapcalc("$output = ($absCal * $input_rast * $esd^2 * $pi) /
($eSun * cos($theta))",
       output=output, absCal=absCalFactor, input_rast=raster,
esd=earthSunDist, pi=math.pi,
       eSun=eSun, theta=theta)

def main():
   imageryPath = r"/path/to/imagery/folder"

   rasterList = grass.list_pairs('rast')

   metadataDict = AssociateMetaDataWithRaster(rasterList)
   rasterVariables = []

   for key in metadataDict:
       grass.message("Extracting variables for %s" % key)
       #Extract variables from xml metadata file given the image name
       #and band number
       variables = ExtractVariablesFromMetadata(os.path.join(imageryPath,
           metadataDict[key]), int(key[-1]))
       #a list of lists containing the raster and variables for TOAR
       #calculation that can be passed to pool.map
       rasterVariables.append([key] + variables)

   pool = multiprocessing.Pool()
   output = pool.map(CalculateTOAR, rasterVariables)
   pool.close()
   pool.join


if __name__ == '__main__':
   #Start multiprocessing of radiometric corrections
   t0 = time.time()/60
   print "Start"
   main()
   print "Finish"
   print (time.time()/60) - t0, "minutes processing time"