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

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!

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])
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")
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)
print('Writing %s for %s'%(ETrast,site))
# Loop through crops