Batch Projection and Masking with GRASS

Submitted by dylan on Wed, 2008-05-07 18:21.

 
Premise
Project and mask 2844 raster data files from the PRISM group, based on a 4km version of CIMIS grids. Good reason to keep a command-line GIS handy.

PRISM Monthly Precip - USAPRISM Monthly Precip - USA

 
Directory Structure

1920-1929
1930-1939
1930-1939
1940-1949
1950-1959
1960-1969
1970-1979
1980-1989
1990-1999

 
Target Grid Geometry (Teal Albers NAD83)

north:      450000
south:      -650000
west:       -400000
east:       600000
nsres:      4000
ewres:      4000
rows:       275
cols:       250
cells:      68750

 
Setup

# create a llwgs72 and teal83 location
#
# longlat WGS72: +proj=longlat +a=6378135 +rf=298.26 +no_defs +towgs84=0.000,0.000,5.000
#
# Teal Albers NAD83: +proj=aea +lat_0=0.0000000000 +lat_1=34.0000000000 +lat_2=40.5000000000 +lon_0=-120.0000000000 +x_0=0.0000000000 +y_0=-4000000.0000000000 +a=6378137 +rf=298.257222101 +no_defs +towgs84=0.000,0.000,0.000 +to_meter=1.0000000000
#
#
# set region resolution to 4km, region extent to sample CIMIS raster, and copy to MASK.

 
Processing

# loop through the list of directories
for dir in 19??-19??
do
echo "Processing $dir files..."

cd $dir

# unzip files (21 seconds per dir)
#
# there may be some bogus files in the lot...
echo "unzipping ..."
for x in *.gz ; do gunzip $x ; done

# llwgs72 location
# load original data (43 seconds per dir):
echo "importing ..."
for x in * ; do r.in.gdal --q -o in=$x out=$x --o 2>/dev/null ; done

# switch to teal83
g.mapset location=teal83 mapset=PERMANENT 2>/dev/null

# teal83 location
# project to teal83 and MASK out everything but CA ( 37 seconds per dir):
echo "projecting and masking ..."
for x in * ; do r.proj --q location=llwgs72  in=$x resolution=4000 out=$x method=nearest --o 2>/dev/null ; done

# save back to disk as ASCII GRID (9 seconds per dir):
echo "exporting to AA GRID ..."
for x in * ; do r.out.arc --q in=$x out=${x}_teal83_ca.asc 2>/dev/null ; done

# archive and compress:
echo "archiving output ..."
tar cfz ../archive_${dir}.tar.gz *.asc && rm *.asc

# clean-up
echo "clean-up ..."
g.mremove  --q -f rast=us* 2>/dev/null

## switch back to llwgs72
g.mapset location=llwgs72 mapset=PERMANENT 2>/dev/null
g.mremove  --q -f rast=us* 2>/dev/null

cd ../

done

PRISM Monthly Precip - CAPRISM Monthly Precip - CA

( categories: )