Applying scalars from a timeseries to raster classes in GRASS GIS (or python)

Question

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!

Answer 1

I was hoping to find a more elegant and less expensive way to do this. I found some examples of folks using pandas with GRASS, but I think they were on OsX or Unix. Getting that to work on my machine wasn't happening anytime soon. Here is my clunky and ugly solution that gives me the right results.

#!/usr/bin/env python


import sys, os, re, fnmatch
import numpy as np
import datetime as dt
from grass.script import parser, run_command, parse_command

def cleanup():
    pass

def main():
    # Path to ETI crop rasters
    rastpath= r"C:\grassdata\idaho_idtm83\IDlandcover\cell"
    # Path to ET depths 
    ETIpath = r"C:\path\to\data"

    # Select only sites with data for the full model run (1985-2015)
    sites = ['BOISE', 'EMMETT','PAYETTE']
    # Raster values in ETI_TV_footprint@IDlandcover 
    sitecodes = [ 1, 2,  3]

    # Get full file path
    files = [os.path.join(ETIpath,'{}_{}'.format(site,'monthly.csv')) for site in sites ]

    # Open files into a dictionary of structured arrays
    d= {} 
    for i,site in enumerate(sites):
        infile = open(files[i])
        flines = infile.readlines()
        infile.close()
        # Extract header for use in numpy array
        hdr = flines[4].split(',')
        # Remove newline
        hdr = [x.strip('\n') for x in hdr] 
        # Put data types and names in tuple  
        dtstr = [(hdr[0], 'datetime64[D]')] + [(x, 'float') for x in hdr[1:]]  
        # Get data
        data = np.genfromtxt( files[i], delimiter=',' , skip_header=5 , dtype= dtstr , names=hdr )
        d[site] = data

    # ----------
    # GRASS STUFF
    # -----------
    # Set region
    run_command("g.region",
                res=30,
                vector="Bound_2010@TV")
    # Mask with TV boundary
    run_command("r.mask",
                overwrite = True,
                vector = "Bound_2010@TV")

And here is where the magic happens

    # 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
            # Extract year from ETI raster
            year = re.search('(\d{4})',ETIrast).group(1) 
            # loop months and make a raster for each month
            for month in range(1,13):
                # loop stations
                ETrast = 'ETdepths_%s_%02d' %( year, month)
                # Initialize zero raster for each month
                expr1 = "tempmap = if( %s@IDlandcover,null())" % (ETIrast)
                run_command("r.mapcalc",
                             overwrite = True,
                             expression = expr1 )
                for i,site in enumerate(d.keys()):  
                    # Boolean array for time indexing the ET dataset
                    idx = d[site]['Date'] == dt.date( int(year),month,1)
                    theseHeaders = d[site].dtype.names
                    print('Writing %s for %s'%(ETrast,site)) 
                    # Loop through crops
                    for crop in theseHeaders[1:]:
                        # TODO Get ET for a crop at site during month during year
                        ET = float(d[site][crop][idx])
                        # TODO creat mapcalc expression iterativelyto create monthly raster
                        expr2 =  "%s = if( %s == %d && ETI_TV_footprint@IDlandcover == %d, %f, tempmap)"  % ( ETrast, ETIrast, int(crop), sitecodes[i], ET) 
                        run_command("r.mapcalc", overwrite=True, expression = expr2)
                        # Copy this map into null map
                        run_command("g.copy", overwrite = True, raster = ETrast +',tempmap')

    return 0

if __name__ == "__main__":
    main()