PerryGeo

Watersheds with lots of roads, buildings, parking lots, rock surfaces, compacted dirt, etc tend to prevent inflitration and cause rapid runoff in response to rainfall. This poses a number of challenges for managing stormwater and water quality. Not surprisingly, the percentage of hydrologically impervious surface in a given watershed is an important factor in many hydrologic models. Using standard aerial photography and GRASS, it’s a relatively simple process to create an impervious surface map using supervised classification.

First find an aerial photo. I grabbed a NAIP image from CASIL but you might want to try using OpenAerialMap. The red, green and blue visible bands are usually sufficient for differentiating between impervious and pervious land use types… For distinguishing different types of vegetation you might want to use a multispectral imagery source with non-visible bands (ie near infrared) but this is usually lower resolution (eg. 30 meter pixels of Landsat) or much more expensive.

Next we jump into GRASS and import our image into a new location:

r.in.gdal -e input=naip.img output=naip location=impervious

Exit and log back into your new location. If you look at the imported rasters, you’ll see three rasters, not one. Each band (R, G and B) gets imported separately.

GRASS 6.3.cvs (impervious):~/>  g.list rast
raster files available in mapset permanent:
naip.1 naip.2 naip.3 

We need to indicate that these rasters form a logical group

i.group group=naip2 subgroup=naip2 input=naip.3@PERMANENT,naip.2@PERMANENT,naip.1@PERMANENT
i.target -c group=naip2

At any time you can list the rasters in a given group/subgroup to confirm.

i.group -l -g group=naip2 subgroup=naip2

Now the real heart of the process. We need to define “training areas” which are polygons around representative land use types. I used QGIS to load the aerial photo and create a new polygon layer with an integer attribute field called vegnum. I digitized a few rocks, paved areas, rooftops and dirt roads to represent the impervious areas to which I assigned vegnum=1. Then I selected some grasslands, forests, lakes and chaparral and assigned 2 as the vegnum. The next step is to load the polygon data into GRASS and rasterize it (in retrospect it would have just been easier to create the grass vector layer from scratch in QGIS to avoid the import step). Note that the vegnum field is specified as the raster value column.

v.in.ogr -o dsn=./training/train1_utm/train1_utm.shp output=train1 layer=train1_utm min_area=0.0001 type=boundary snap=-1
v.to.rast input=train1 output=train1 use=attr column=vegnum type=point,line,area layer=1 value=1 rows=4096

Next we use i.gensig to generate a spectral signature (the statistical profile; mean and covariance matrix of the input pixels) for the training areas.

i.gensig trainingmap=train1 group=naip2 subgroup=naip2 signaturefile=naip2_train1.sig

Now that we have a signature of impervious vs. non-impervious surfaces, we can use the maximum likelihood method to classify each pixel into the highest probability category.

i.maxlik group=naip2 subgroup=naip2 sigfile=naip2_train1.sig class=imperv

You might notice a slight speckled, noisy appearance due to things like shadows, reflections or imperfect training areas. Usually these small 1-pixel deviations are not interesting enough to keep so we can smooth out the image taking the mode (most comon) cell in a 3×3 window.

r.neighbors input=imperv output=imperv_mode method=mode size=3 

And here are the results… calculating imperviousness will most likely be an iterative process so be prepared to evaluate the output, tweak the training areas and rerun the process a few times. Once you’re happy with the results, you can use zonal statistics with a tool like starspan to find the percent imperviousness of your watersheds or other regions.

It just keeps getting easier and easier to get a fully functional open source GIS workstation up and running thanks to Ubuntu. The following instructions will take your vanilla installation of Ubuntu 7.10 and add the following top-notch desktop GIS applications:

• Postgresql/PostGIS : a relational database with vector spatial data handling
• GRASS : A full blown GIS analysis toolset
• Quantum GIS: A user-friendly graphical GIS application
• GDAL, Proj, Geos : Libraries and utilities for processing spatial data
• Mapserver : web mapping program and utilites
• Python bindings for QGIS, mapserver and GDAL
• GPSBabel : for converting between various GPS formats
• R : a high-end statistics package with spatial capabilities
• GMT : the Generic Mapping Tools for automated high-quality map output

While this is not a comprehensive list of open source GIS software, these packages cover most of my needs. If you want to live on the bleeding edge and have to have the absolute latest versions, you’ll be better off installing these from source. But for those of us that want a stable and highly functional GIS workstation with minimal fuss, this is the way to go:

1. Go to System > Administration > Software Sources and make sure the universe and multiverse repositories are turned on. Close the window and the list of available software packages will be refreshed.
2. Open up a terminal (ie the command line) via Applications > Accessories > Terminal and type the following:
   

   sudo apt-get -y install qgis grass qgis-plugin-grass mapserver-bin gdal-bin cgi-mapserver \
   python-qt4 python-sip4 python-gdal python-mapscript gmt gmt-coastline-data \
   r-recommended gpsbabel shapelib libgdal1-1.4.0-grass

   The sudo part indicates that the command will be run as the administrator user, apt-get -y install is the command telling it to install the list of packages and answer yes to any questions that pop up.

3. There is one package that is worth upgrading to the latest and greatest - Quantum GIS. The latest version (0.9) is due out very shortly and has the ability to write plugins using the python programming language. A big plus!

   Download the latest build from http://qgis.org/uploadfiles/testbuilds/qgis0.9.0.debs_ubuntu_gutsy.tar.gz and extract it ( right-click > Extract Here ). In the directory you’ll see 4 .deb files, only 3 of which you’ll need unless you plan on doing any development work.

   Double click libqgis1_0.9.0_i386.deb and you’ll get a message saying an older version is available from directly from ubuntu. We already know this so just close and ignore it. Click Install Package and wait for it to complete then close out.

   Repeat for qgis_0.9.0_i386.deb and qgis-plugin-grass_0.9.0_i386.deb (in that order).

And there we have it, about 15 minutes depending on your internet speed and you’ve installed a high-end GIS workstation built completely on free and open source software.

LIDAR data is certainly a hot technology these days. LIght Detection And Ranging data can be used to create extremely detailed terrain models but there are lots of barriers to using LIDAR data effectively. USGS Center for LIDAR Information Coordination and Knowledge was put in place to “facilitate data access, user coordination and education of lidar remote sensing for scientific needs“.

Beyond the sheer size of the datasets and the knowledge and hardware required to process them, software is a big issue. In the realm of open-source GIS tools, there are many applications (GRASS being the most prominent) for dealing with elevation point data and processing it into more meaningful products such as elevation DEMs and contours.

Usually the data comes as simple ASCII text files and the x,y and z values are easily extracted from such a file. But take a look at the USGS data distribution site and you’ll notice some of the datasets are distributed as LAS binary files. It makes sense to store such massive datasets in binary so I started looking for some LAS conversion tools.So after some searching, I found a bunch of proprietary products for working with LAS but no open source tools. Luckily, the format is well documented thanks to the efforts by the ASPRS to make it an open specification.

So dusting off my notes about parsing binary files in python, I set out to create a python module for extracting LIDAR data from LAS files. The LAS format contains a header which needs to be parsed first in order to read the point cloud. Once you have the header info, you can scan your way through the dataset to pick out the x,y,z values.

Here’s an example of the python interface that will read the first 10,000 points into a 2D shapefile with the elevation as a attribute in the dbf:

import pylas
infile = 'sanand000001.las'
outfile = 'lidar.shp'
header = pylas.parseHeader(infile)
pylas.createShp(outfile, header, numpts=10000, rand=False)

The issue I struggled with is the sheer size of these datasets. A USGS quarter quad can contain 10 million points which is an excessive number of points to create, say, a 10 meter DEM over such a small area. Clearly there was a need to extract a subset of this dataset but just taking the points sequentially gives you a subset of the total area. So, by default, pylas randomly scans the data to pull the number of specified points so that the point cloud could cover the entire area (at a much lower point density). Without numpts specified, it will randomly select 1/2000th of the total number.

So the simplified interface to make a more manageable lidar shapefile would be:

header = pylas.parseHeader(infile)
pylas.createShp(outfile, header)

Once the shapefile is created, you can bring it into GRASS to do the processing to generate DEMs, contours and other derived elevation products:

v.in.ogr dsn=lidar.shp layer=lidar output=lidar
g.region vect=lidar
g.region res=10
v.surf.rst input=lidar elev=lidar_dem zcolumn=elev

# Launch the interactive 3D viewer
nviz lidar_dem

Of course the method I just described is very simplistic and does not even come close to utilizing the full potential of the LIDAR point cloud, but it’s a start.

The pylas.py module can be downloaded here. The code has worked for me on the few datasets I’ve tested it with but it should certainly be considered a rough-cut, alpha product. There is much room for improvement and, of course, if you have any suggestions or contributions, please get in touch.

NOAA Electronic Navigational Charts (ENC) contain (among many other things) depth soundings that can be processed into raster bathymetry grids. The ENC files are available as a huge torrent from geotorrent.org (http://geotorrent.org/details.php?id=58).

Download this torrent and check readme.txt to find the chart of interest:

Port Hueneme to Santa Barbara|5|2005-10-03|2005-10-03|US5CA65M

First check out the gdal documentation for s57 files at http://www.gdal.org/ogr/drv_s57.html.

Change to the US5CA65M directory and you’ll see a .000 file (and maybe .001, .002 etc). Run ogrinfo on the .000 file and you’ll see ~ 61 layers, one of which (”SOUNDG”) represents the soundings. Let’s start by examining the soundings layer:

ogrinfo -summary US5CA65M.000 SOUNDG

We see that there are 43 “features” but since the features are multipoints, there are actually thousands of soundings. The multipoints are 3D so If we convert to a shapefile with ogr2ogr’s default settings we loose the 3rd dimension. To solve this, we need to append “25D” to the layer type. Furthermore, the multipoint geometry confuses some applications so we want to split it into a layer with simple 3D point geometries. Luckily there is a SPLIT_MULITPOINT option that must be specified as an environment variable:

export OGR_S57_OPTIONS=”RETURN_PRIMITIVES=ON,RETURN_LINKAGES=ON,LNAM_REFS=ON,SPLIT_MULTIPOINT=ON,ADD_SOUNDG_DEPTH=ON”
ogr2ogr -nlt POINT25d test3.shp US5CA65M.000 SOUNDG

Now we get ~ 3000 3D points with the depth added as an attribute for good measure.

Now bring these into grass and create a raster:

v.in.ogr -zo dsn=test3.shp output=soundg layer=test3
v.info soundg
g.region vect=soundg nsres=0.001 ewres=0.001
v.surf.rst input=soundg elev=bathy layer=0
r.info bathy

since depths actually show up as positive elevations, we want to multiply the grid by -1

r.mapcalc sb_bathy=bathy*-1

And of course we want to make some nice shaded relief and contour maps for viewing with QGIS:

r.shaded.relief map=sb_bathy shadedmap=sb_shade altitude=45 azimuth=315
r.contour input=sb_bathy output=sb_contour step=5
qgis &

From the screenshot, we see the pits and spikes from potential outliers so we might want to go back and adjust the tension and smoothing on the raster creation (the v.surf.rst command).