Further adventures with AHPS and HRAP

So, once I figured out how to deal with the oddball projection, I knew in principle how to deal with precipitation shapefiles provided by NOAA. However, going through each file one-by-one would be a pain. I need to automate this, or at least partially automate it. QGIS comes with a graphical modeler that allows limited automation of tasks, but I quickly found a problem with this–I need to do some reprojection of layers from one custom CRS to another, and the “reproject layer” function available in the graphical modeler does not work with custom CRS.

Fortunately, poking around online I found a simpler solution–use shell scripts and the command line in Terminal. Fortunately, the reprojection tool used in QGIS (ogr2ogr) handles custom CRS just fine–you can just specify the projection in proj4 format and away it goes. So, now all I needed to do is figure out how to get syntax for the parameters and whatnot to work. Luckily, QGIS provides a transcript of what commands it is sending when you run a tool from one of the menus, so I can (mostly) just deal with QGIS and copy-paste from there to an sh file. I ended up with producing a set of four folders (one for the raw shapefiles from AHPS and one for each stage in the process) and three shell scripts, each of which does some needed modification of the layers and puts the output in the next folder. Here are the sh files, first, for reprojecting the AHPS files from the HRAP geographic coordinate system to the HRAP projected coordinate system:

for f in *.shp
do
  echo "Processing $f"
/Library/Frameworks/GDAL.framework/Versions/1.11/Programs/ogr2ogr 
-s_srs "+proj=longlat +a=6371200 +b=6371200 +no_defs" -t_srs 
"+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-105 +k=1
+x_0=0 +y_0=0 +a=6371200 +b=6371200 +to_meter=4762.5 +no_defs" 
../reprojected/$f $f
done

And then rasterization:

for f in *.shp
do
  echo "Processing $f"
/Library/Frameworks/GDAL.framework/Versions/1.11/Programs/gdal_rasterize 
-a Globvalue -tr 1.0 1.0 $f ../rasterized/$f.tif
done

Last, clipping the rasters to a smaller extent so that the collection of files is not unnecessarily large:

for f in *.tif
do
  echo "Processing $f"
/Library/Frameworks/GDAL.framework/Versions/1.11/Programs/gdal_translate 
-of GTiff -projwin -112.0 -1327.0 17.0 -1414.0 $f ../output/$f
done

That all works nicely and allows me to process large numbers of files quickly and efficiently. I’m sure that if I really understood sh (rather than copying the work of others and modifying things without really knowing what I’m doing) I could probably just have an “input” and “output” directory and a single sh file. But this is still pretty efficient and didn’t require too much farting around.

Leave a Reply