Extracting the basin boundary from subbasins

Assume that you have a number of adjacent subbasins (e.g. the USGS HUCs) but your study area is a superset of those. All your subbasins are contained within one vector map, or if they aren’t you have patched them together and selected the subbasins you need. Now we want to get the boundary for plotting a nice map of the study area or easily getting a model grid for example. In GRASS GIS, we use

v.clean input=map output=mapc tool=snap,break,rmdupl thresh=.01
and
v.dissolve input=mapc output=mapb col=reg
where map is your original vector map, reg is just the common column that you want your inner boundaries dissolved by. Now I just need to figure out how the threshold value should change in the v.clean command.

Generating random fields using Rpy2

Following on the previous post, here’s one way to generate random fields using Rpy2 and the R RandomFields package. Let’s assume your data are in the lat, lon, data vectors, first we import modules, set some options and fit a variogram (check the RandomFields documentation for details)

import rpy2.robjects as robjects
import rpy2.robjects.numpy2ri
R=robjects.r
R.library('RandomFields')
TF={'table.format':True,'lsq.methods':'self','mle.methods':'NULL'}
estparam=robjects.FloatVector((mean(data),var(data),None,None,None))
F=R.fitvario(lon,lat,data=datam,model='whittlematern',param=estparam,**TF)

The number of parameters depends on the model selection, having the format (mean,variance,nugget,scale,…), while each parameter that needs to be estimated is set to None when calling fitvario. The option table.format sets the output to a matrix, so now we can get the parameters (first index is column, second is row), do a

print(F)

to check what the available values are, e.g.

nugget=F[2][3]
scale=F[2][2]
kappa=F[2][1]

Then we can simulated random field realizations

estparam=robjects.FloatVector((mean(data),var(data),nugget,scale,kappa))
newdata=R.GaussRF(lon,lat,grid=False,model='whittlematern',param=estparam,n=num)

where num is the number of realizations.

Inverse distance weighting interpolation in Rpy2

One thing that’s nice about Python is that you can easily use it to glue lots of interesting software together. Take Rpy for example: calling R functions from within python. Let’s say we want to perform IDW interpolation from spatial data. I will assume here that lon,lat are vectors with the locations of the existing data, and vlon,vlat the locations of the data to be interpolated, while data is a vector that contains your (wait for it) data.
First, let’s import the modules

import rpy2.robjects as robjects
import rpy2.robjects.numpy2ri
R=robjects.r
R.library('gstat')

Let’s create the data frames, the first one containing your measurements

d={'x':lon,'y':lat,'data':data}
D=R['data.frame'](**d)

and then create a data frame containing the locations where the data will be interpolated

dn={'x':vlon,'y':vlat}
DN=R['data.frame'](**dn)

and let’s run the IDW

ld=R.idw(robjects.RFormula('data~1'),robjects.RFormula('~x+y'),D,DN)

and then extract the interpolated data

newdata=array(ld.subset('var1.pred')[0])


Importing SNODAS data into GRASS GIS

From the SNODAS website: “…is a modeling and data assimilation system developed by NOHRSC to provide the best possible estimates of snow cover and associated variables to support hydrologic modeling and analysis. The aim of SNODAS is to provide a physically consistent framework to integrate snow data from satellite and airborne platforms, and ground stations with model estimates of snow cover (Carroll et al. 2001)”. So, how can we import it into the GRASS GIS?

First download the data, and untar them (once for each month, and once for each day), and you should get pairs of “.dat” and “.Hdr” files. The data files are stored in flat 16-bit binary format, so assuming that “snowdas_in.dat” is the name of the input file, at the GRASS prompt:

r.in.bin -bs bytes=2 rows=3351 cols=6935 north=52.874583333332339 south=24.949583333333454 east=-66.942083333334011 west=-124.733749999998366 anull=-9999 input=snowdas_input.dat output=snowdas

The units are described in the SNODAS documentation (section 2 table), so use r.mapcalc to convert the data.

Extracting SNOTEL SWE measurements

Here’s a simple Python script to extract SNOTEL SWE measurements from the inconveniently formatted files that can be downloaded at NRCS.


#!/usr/bin/python
# Extracts SNOTEL SWE observations from NCDC-format to standard ASCII-column format
# import modules
import sys
import string
import datetime
# get and open filename
fid=open(sys.argv[1])
# read file contents
ln=fid.readline()
while (ln!='' and string.find(ln,'Snow')>-1):
yr=int(string.split(ln)[1])+1899
if (yr<1910): yr+=100
sdate=datetime.date(yr,10,1)
dy=0
fid.readline()
fid.readline()
fid.readline()
fid.readline()
fid.readline()
fid.readline()
lns=[]
for t in range(0,31):
lns.append(fid.readline())
for m in range(0,12):
for ln in lns:
s=string.strip(ln[(4+m*6):(10+m*6)])
if (s!='---'):
date=sdate+datetime.timedelta(dy)
dy+=1
if (s):
print('%04d %02d %02d %.2f' % (date.year,date.month,date.day,float(s)*25.4))
else:
print('%04d %02d %02d %.2f' % (date.year,date.month,date.day,-99.))
fid.readline()
fid.readline()
fid.readline()
fid.readline()
fid.readline()
ln=fid.readline() ## WHILE LOOP ENDS HERE


This will output a nice ASCII file with YEAR MONTH DAY SWE in its columns (SWE is in mm).

Plotting maps in Python

Here’s a simple script making use of Matplotlib and their mapping toolkit Basemap to make a quick map of your data, assuming that you have them as a vector along with their latitude, longitude as points. Obviously you can change the cylindrical projection used here (take a look at the Basemap documentation for available options).


# Map plotting function using Basemap (Cylindrical projection)
# Input is LON, LAT, DATA, (MINLON,MAXLON,MINLAT,MAXLAT,RES)
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as Pyplot
import numpy as Numpy
def plotmap(lon,lat,data,gd):
# define domain from gd tuple
ll_lon=gd[0]
ll_lat=gd[2]
ur_lon=gd[1]
ur_lat=gd[3]
# set grid definition
res=gd[4]
cols=int((ur_lon-ll_lon)/res)+1
rows=int((ur_lat-ll_lat)/res)+1
lon0=ll_lon
lat0=ll_lat
# assign data
grid=Numpy.ma.masked_all((rows,cols))
for c in range(0,len(data)):
i=int((lat[c]-lat0)/res)
j=int((lon[c]-lon0)/res)
grid[i,j]=data[c]
# create new figure
fig=Pyplot.figure()
ax=fig.add_subplot(1,1,1)
# put background map and plot data
m=Basemap(projection='cyl',llcrnrlon=ll_lon,llcrnrlat=ll_lat,urcrnrlon=ur_lon,urcrnrlat=ur_lat,ax=ax)
im=m.imshow(grid,interpolation='nearest')
fig.colorbar(im)
m.drawstates()
m.drawrivers()
Pyplot.show()

I adapted this from someone off some mailing list, but can’t remember where, so thanks to that guy/gal.

Reprojecting MODIS snow cover data to a grid

MODIS snow cover extent satellite data (products MOD10A and MYD10A) are provided in the Integerized Sinusoidal (ISIN) projection at 500 m spatial resolution (among others). I was looking for a quick way to reproject the data to a standard geographical grid, such as the one used by the Variable Infiltration Capacity (VIC) model.
First, download the MODIS Reprojection Tool (MRT) and install it according to instructions (the binary version seems to work fine). Now order the data from WIST, selecting the dates and the spatial region. This will give you one or more tiles (1200×1200 km) of data in HDF-EOS format.
Now having multiple tiles, we have to mosaic them.

mrtmosaic -i file_list -o mosaiced_file.hdf

where file_list contains the name of the files to be mosaiced. Then resample using a parameter file created from the ModisTool GUI for example if you don’t want to crate it manually.

resample -i mosaiced_file.hdf -o final_file.hdf -t GEO -p parameter_file

We now have our spatially subsetted (sic) HDF file which we can process with a simple python script

#!/usr/bin/python
# Read MODIS SCE data and aggregate to a 1/8th degree VIC grid
import numpy
import os
import string
import re
from pyhdf.SD import SD
# get VIC grid
x=numpy.loadtxt(VIC SOILFILE)
vlat=x[:,2]
vlon=x[:,3]
# read MYD files
t=1
mod=SD('../MODIS/MYD10A2.A2005%03d.hdf' % t)
sc=mod.select('Eight_Day_Snow_Cover')
msc=mod.select('Maximum_Snow_Extent')
at=mod.attributes()
s=at['OldArchiveMetadata.0']
S=string.split(s,"NORTHBOUNDINGCOORDINATE")
maxlat=float(re.findall("d+.d+",S[1])[0])
S=string.split(s,"SOUTHBOUNDINGCOORDINATE")
minlat=float(re.findall("d+.d+",S[1])[0])
S=string.split(s,"EASTBOUNDINGCOORDINATE")
maxlon=float(re.findall("-*d+.d+",S[1])[0])
S=string.split(s,"WESTBOUNDINGCOORDINATE")
minlon=float(re.findall("-*d+.d+",S[1])[0])
xdim=sc.dimensions()['XDim:MOD_Grid_Snow_500m']
ydim=sc.dimensions()['YDim:MOD_Grid_Snow_500m']
lat=numpy.zeros((ydim,xdim))
lon=numpy.zeros((ydim,xdim))
for i in range(0,ydim):
for j in range(0,xdim):
lat[i,j]=maxlat-(maxlat-minlat+1.0)/float(ydim)*i
lon[i,j]=minlon+(maxlon-minlon+1.0)/float(xdim)*j
lat=numpy.reshape(lat,xdim*ydim)
lon=numpy.reshape(lon,xdim*ydim)
sce=numpy.reshape(sc.get(),xdim*ydim)
msce=numpy.reshape(msc.get(),xdim*ydim) for c in range(0,len(vlat)): idx=numpy.where(numpy.where(numpy.where(numpy.where(lat(vlat[c]-0.125/2.),lon,numpy.NaN)>(vlon[c]-0.125/2.),lon,numpy.NaN)<(vlon[c]+0.125/2.))[0]
# find snow
snow=float(len(numpy.where(msce[idx]==200)[0]))/float(len(idx))
# find land
land=float(len(numpy.where(msce[idx]==25)[0]))/float(len(idx))
# find cloud
cloud=float(len(numpy.where(msce[idx]==50)[0]))/float(len(idx))
print('%.4f %.4f %.3f %.3f %.3f' % (vlon[c],vlat[c],snow,land,cloud))

That should output your regridded MODIS snow cover data in a GMT-ready format. Of course, the regridding method is not entirely accurate, since the aggregate grid cell values are not the area-weighted averages; instead, they are computed from the 500-m MODIS pixels that are contained within each coarser grid cell.

Hello world!

Hmmmm, what the title says. Mostly a little corner to document some of the stuff I do with my research (that might be useful to somebody else as well), and rant about my favorite music, films and other non-serious stuff…