I have some yearly crop data layer rasters to which I would like to apply monthly crop-specific evapotranspiration depths from a text file. I could go month by month and do some raster math, but that seems like the clunky way to go about it. There are also several weather stations with different ET values for each crop, so my clunky expression would be something like
r.mapcalc expression = "ETdepths2015 = if( CDL_2015@IDlandcover == 1 & stationfootprint@IDlandcover == 12, 0.3) + if( CDL_2015@IDlandcover == 4 & stationfootprint@IDlandcover == 12, 0.4) + ..."
Problem is there are 50 unique crop class, 12 months in the year, 10 years, and 9 stations, so that would be about 54000 lines of mapcalc
I'm new to scripting in GRASS and I would like to use numpy to do the calculations without having to export the rasters, but I haven't gotten to far in that regard. I'll keep working on it, though if anyone has any hints, I'd be happy to hear them.
Update and some ideas
I was hoping use grass.script.array.array() to convert the rasters to numpy arrays and match crop codes, though I couldn't get that to work.
Here is a sketch of how I think I could still use r.mapcalc by building an expression iteratively.
# Loop ETI years for root,dirs,fnames in os.walk(rastpath): for fname in fnmatch.filter(fnames,'ETI_Idaho*'): ETIrast = fname # String of raster map # TODO get year from raster map year = 2015 # loop weather stations for site in d.keys(): # loop through months and make a raster for each month for month in range(1,13): # Boolean array for grabbing specific row idx = dates == dt.date(year,month,1) # Loop through for crop in hdr[1:]: # Get ET for a crop at site during month during year ET = d[site][crop][idx] # Mapcalc expression to create monthly raster expr = '"ETdepths_%s_%s" = if( ETI_TV_voronoi@IDlandcover == %s & ETIrast == %s, %s ) + ' % (year, month,site,crop, ET)
And....I'm stuck on how to do that iterative string building!