Raster Operations
Geologic-Scale Erosion
100 iterations of mass removal based on preferential flow of water as calculated by r.topidx in GRASS. Notice how landform is cutdown most in stream channels, least at the ridges. The basins between ridges appear to "fill" with sediment near the 50th iteration as the entire landform is lowered to sea level (0m). Absolute change in elevation is visble in the elevation profiles below. Example GRASS commands below. Here is a link to a movie, containing all 100 iterations.
Note that this is *not* a direct measure of sediment removal, mearly one means of visualizing the qualitative effects of sediment removal. Actual sediment flow values could possibly be calculated with r.terraflow and some expert knowledge of the area. In addition, sediment transport models such as RUSLE, USPED, and SIMWE will provide a quantitative erosion AND deposition estimate. Details here:
- Comparision of the three methods by H. Mitasova
- Explanation of SIMWE by H. Mitasova
- Explanation of USPED by H. Mitasova
- Additional erosion modeling information by H. Mitasova
- http://skagit.meas.ncsu.edu/~chris/terradyn/lw.html
#
for x in `seq 0 1 100`
do r.mapcalc "e = if( float(elev_meters - float($x*ti)) > 0.0 ,float(elev_meters - float($x*ti) ), 0.0 )"
r.shaded.relief map=e shade=s1
d.rast s1
d.barscale at=0,0
d.out.png out=e_$x res=1
d.erase
done
convert `ls *.png | sort -n -k2 -t"_"` movie.mpg
#make an elevation profile
r.profile in=e100 profile=663451.38153808,4038227.40942676,663719.16510544,4041226.58538121 > e100.profile
Mapping Wifi Networks with Kismet, GDAL, and GRASS
Some simple examples of how to extend data collected from kismet along with gpsd. The prgram gpsmap, bundled with kismet, can create excellent maps of wifi networks, albeit not in a georeferenced format. With some ideas from Matt Perry, a quick read through some python docs, and a little testing I put together a simple python script (kismet2shp) which will convert Kismet XML logfiles into a shapefile (see attached file at bottom of page for the source). Conversion from kismet logfile to KML was demonstrated in the above article by Matt. Simple usage of kismet2shp is summarized:
./kismet2shp -r -s --file=Kismet-Sep-03-2006-1.xml --outfile=range.shp
#
# network centers (points)
./kismet2shp -s --file=Kismet-Sep-03-2006-1.xml --outfile=points.shp
Output can either be in the default latlong WGS84 GCS, or in some user defined projected coordinate system (see source code). Example images from gpsmap and from QGIS demonstrate similarities. The important difference is that by converting Kismet data into a shapfile, it is possible to bring this information into a GIS. A simple example below demonstrates how one can compute an estimated signal power index in GRASS, similar to the interpolated power option in gpsmap. The power density index is calculated as the summation of [ 1/(range^2) ] for each network center.
gpsmap range: network range as visualized with gpsmap
gpsmap power: interpolated power density from gpsmap
qgis range: network range shapefiles visualized in QGIS
GRASS power density: simple index of wifi signal power computed in GRASS
Make simple maps with gpsmap
gpsmap -S 2 -e -n 1 -t -l ssid -k -p -Q 20 -q 2 -P 30 -z 2 -d 1000,1000 \
-o map1.png Kismet-Sep-03-2006-1.gps
#range map
gpsmap -e -n 1 -t -l ssid -k -r -R 15 -q 2 -d 1000,1000 \
-o map1.png Kismet-Sep-03-2006-1.gps
Estimate signal power density in GRASS
#import data from kismet2shp:
v.in.ogr dsn=. layer=points out=wp
#setup a resonable resolution
g.region vect=wp res=10 -a
#init a constant value surface for distance calculation
r.mapcalc "z = 1"
#init power level map
r.mapcalc "a = 0.000"
#
#create a list of networks
cats=`v.db.select wp -c colum=cat`
#loop thru list of networks
for i in $cats
do
v.extract in=wp out=wp_i list=$i --o > /dev/null
#get this network's estimated range
range_i=`v.db.select wp_i -c colum=range`
#compute distance map from network center to estimated range
r.cost -k in=z out=c_i start_points=wp_i max_cost=$range_i --o > /dev/null
#convert distance from net center to power index
r.mapcalc "c_i_f = if( isnull(c_i) != 1, 1/(float(c_i)^2) , null())" > /dev/null
#add to previous power estimation
r.mapcalc "a = if( isnull(c_i_f) != 1, a + c_i_f, a)" > /dev/null
done
#remove the weakest power levels
r.mapcalc "a = if( a < 0.5, null(), a)" > /dev/null
#fix colors
r.colors map=a color=gyr
A mobile wifi mapping setup might include:
Planned work:
- Modifications to kismet-earth to improve merging of multiple .gps files
- enhanced wifi network center estimation, and power interpolation with GIS tools
Simple Comparision of Two Least-Cost Path Approaches
Least Cost Path Comparision Map r.walk in blue, r.cost in orange
Premise:
Simple comparison between the least-cost path as computed by r.walk vs. r.cost. Default parameters, unity travel friction, and λ = 1 were used with r.walk. Slope percent was used as the accumulated cost for r.cost. Test region similar to that used in a previous navigation example.
Results:
- r.walk
- Shorter total path (6.03 km)
- Tendency to stay mid-slope, and cross large drainage at the gentlest possible position
- r.cost
- Longer total path (7.12 km)
- Tendency to follow drainages upstream-- not always possible in steep terrain
Code:
Prep work see attached files at bottom of page to re-create
# unity friction for r.walk
r.mapcalc "friction = 1.0"
# use slope percent as the friction for r.cost
r.slope.aspect elev=elev slope=slope format=percent
Compute the least-cost path via two methods
r.walk -k elev=elev friction=friction out=walk.cost start_points=start stop_points=end lambda=1
# Peak cost value: 14030.894105
r.cost -k in=slope out=cost start_points=start stop_points=end
# Peak cost value: 11557.934493
# compute shortest path from start to end points
r.drain in=walk.cost out=walk.drain vector_points=end
r.drain in=cost out=cost.drain vector_points=end
Vectorize least-cost paths for further analysis
r.to.vect -s in=cost.drain out=path_2
# add DB column for length calculation
v.db.addcol path_1 columns="len double"
v.db.addcol path_2 columns="len double"
# how long are the two paths?
v.to.db path_1 option=length column=len units=me
# 6.03 km
v.to.db path_2 option=length column=len units=me
# 7.12 km
Simple map
d.vect start icon=basic/circle fcol=green size=15
d.vect end icon=basic/circle fcol=red size=15
d.vect path_1 width=2 col=blue
d.vect path_2 width=2 col=orange
Sample the elevation raster at a regular interval
# this will create two tables for each
# path_1_pts_1 <-- original atts
# path_1_pts_2 <-- original cat + distance along the original line
v.to.points in=path_1 out=path_1_pts -i dmax=100
v.to.points in=path_2 out=path_2_pts -i dmax=100
# add column for elevation values
v.db.addcol map=path_1_pts layer=2 columns="elev double"
v.db.addcol map=path_2_pts layer=2 columns="elev double"
# upload elev values for profile graph
v.what.rast vector=path_1_pts raster=elev layer=2 column=elev
v.what.rast vector=path_2_pts raster=elev layer=2 column=elev
# dump profiles to CSV files
db.select path_1_pts_2 fs="," > p1.csv
db.select path_2_pts_2 fs="," > p2.csv
Least Cost Path Comparision Profile
Plot an elevation profile for each path
R
## load data from previous steps
p1 <- read.csv('p1.csv')
p2 <- read.csv('p2.csv')
## can't guarantee that the two lines will be in the same direction (a r.to.vect thing?)
## p1 is correct, p2 is reversed
## make a nicer plot
## library(Cairo)
## Cairo(type='png', bg='white', width=800, height=400)
plot(elev ~ rev(along), data=p2, type='l', col='orange', lwd=2, xlab='Travel Distance (meters)', ylab='Elevation Change (meters)')
lines(elev ~ along, data=p1, col='RoyalBlue', lwd=2)
legend(0, 3000, legend=c('r.walk', 'r.cost'), col=c('RoyalBlue', 'orange'), lwd=2)
Visual Comparison of 2 Raster Images
The visual comparison of raster data can be complicated by images that have non-NULL values at each grid cell. Transparency can used to "see" both images at once, but the results can be misleading. An alternative method suggested in the GRASS Book (Neteler and Mitasova, 2003), uses a couple simple raster algerbra tricks to make a "checkerboard" image composite of both input raster files. Figure 1 illustrates a simple comparison between an IKONOS false infra-red image, and a classified vegetation image via this checkerboard approach.
Example Code
r.mapcalc " v=if( (sin(5*row()) + sin(5*col())) > 0.0, null(), veg_11 )"
r.mapcalc " i=if( (sin(5*row()) + sin(5*col())) > 0.0, ikonos.rgb, null() )"
#assign the original colors to the new files:
r.colors map=v rast=veg_11
r.colors map=i rast=ikonos.rgb
#patch them together with r.patch, or view with:
d.rast v
d.rast -o i
References:
Neteler, M. & Mitasova, H. Open Source GIS A GRASS GIS Approach. Kluwer Academic Publishers, 2003.
Watershed Basin Delineation
Multiple approaches to the automatic delineation of watershed boundaries, along with the processing associated with the CalWater 2.2.1 dataset. Note that the output from each algorithm is slightly different, and that CalWater does not always represent physiographicaly correct boundaries. The CalWater 2.2.1 original dataset and converted to shapefile are available for download. See attached DEM clip for examples presented below. Matt Perry has written up a nice tutorial on identifying impervious surfaces within the context of watershed delineation.
Watershed basin delineation by the following approaches:
- ArcInfo basin() GRID function
- GRASS r.watershed command
- GRASS r.terraflow command
multiscale approach suitible for small to medium input grids. raster output.
minimum basin size must be defined suitible for small to medium input grids. raster output.
multiscale approach suitible for massive input grids, fast. raster output.
Notes
- none of these algorithms can deal with hierarchial watershed decomposition
- removal of small vector features does not ensure perfect basin delineation
- raw basin output from these algorithms are best further processed some ideas here
watershed 1: 3D comparison of the four methods. calwater (red), basin (black), r.terraflow (white).
watershed 2: 3D comparison of the four methods. calwater (red), basin (black), r.terraflow (white).
Pre-processing in GRASS note that these steps can be done in any GIS
# set the region and cut out a sample section of the DEM
g.region res=90 -ap
r.mapcalc "a = ca_90m_dem_meters"
#generate a shaded relief map for viz:
r.shaded.relief map=a shade=a.shade
ArcInfo basin()
# note that 'a_grid' was exported from GRASS via: r.out.arc in=a out=a_grid
asciigrid a_grid a float
grid
basins = basin(flowdirection(a))
GRASS r.watershed
# note that a threshold of 200 sets the minimum basin size [update this]
r.watershed --o elev=a threshold=200 basin=w.basins
GRASS r.terraflow
r.terraflow elev=a filled=a.filled direction=a.dir swatershed=a.swater accumulation=a.acc tci=a.tci
#set cell values of 0 to NULL, these are regions where calculation was not possible
r.null map=a.swater setnull=0
Post-processing of generated basins in GRASS note that this isn't a true merging of sub-basins, just a rough idea of what it might look like!
# r.terraflow
r.to.vect -s in=a.swater out=a_watersheds feature=area
# ArcInfo basin()
r.to.vect -s in=arc.basins out=arc_basins feature=area
# remove basins smaller than approx. 200ac.
# r.terraflow
v.clean --o in=a_watersheds out=b_watersheds tool=rmarea thresh=809371
# ArcInfo basin()
v.clean --o in=arc_basins out=arc_basins_clean tool=rmarea thresh=809371
CalWater note that this comes as an e00 file, and must be converted.
# note that this is an e00 file and needs to be converted
# converted with ArcInfo
import AUTO calw221.e00 calwater
# converted with AVCEimport http://avce00.maptools.org/avce00/index.html
~/src/avce00-1.3.0/avcimport calw221.e00 calwater_avc
# ArcMap simple conversion to shapefile and re-projection possible via standard tools
# converted to shapefile polygon
# conversion to shapefile with attributes intact not as simple...
# projection to our spatial ref system via ogr2ogr
# note that this shapfile was created by ArcMap
ogr2ogr -t_srs '+proj=aea +x_0=0.0 +y_0=0 +lon_0=-96 +lat_0=40.0 +lat_1=20 +lat_2=60.0 +datum=NAD83' calwater_aea.shp calwater\ polygon.shp
Visual comparison in GRASS
v.in.region out=temp_region
v.overlay ainput=calwater binput=temp_region out=calwater_clip operator=and
g.remove vect=temp_region
## compare the methods:
d.his i=a.shade h=w.basins # r.watershed
d.vect calwater_clip type=boundary col=red width=2 # Calwater
d.vect arc_basins_clean type=boundary col=black width=1 # ArcInfo basin()
d.vect b_watersheds type=boundary col=white width=1 # r.terraflow
Aggregation of CalWater Data
#aggregation of CalWater data: see metadata:
Relevant columms:
43 HRC 2 2 C - Hydrologic Region Code (DWR)
45 HBPA 2 2 C - Hydrologic Basin Planning Area (RWQCB)
47 RBU 5 5 I - Concatenates HR,RB,HU into a single integer
52 RBUA 6 6 I - Concatenates HR,RB,HU,HA
58 RBUAS 7 7 I - Concatenates HR,RB,HU,HA,HSA
65 RBUASP 9 9 I - Concatenates HR,RB,HU,HA,HSA,SPWS
74 RBUASPW 11 11 I - Concatenates HR,RB,HU,HA,HSA,SPWS,PWS
85 HR 2 2 I - Hydrologic Region (1->10)(DWR)
87 RB 1 1 I - Regional Water Qual. Cont. Board (1->9)(RWQCB)
88 HU 2 2 I - Hydrologic Unit (00->~80)(SWRCB)
90 HA 1 1 I - Hydrologic Area (0->9)(SWRCB)
91 HSA 1 1 I - Hydrologic Sub-Area (0->9)(SWRCB)
92 SPWS 2 2 I - Super Planning Watershed (00->~30)(CDF)
94 PWS 2 2 I - Planning Watershed (00->~13)(CDF)
# polygons can be dissolved via ArcMap / PostGIS / v.to.rast use=attr -> r.to.vect
regarding watershed decomposition and limitation C/O Helena Mitasova
Dylan, r.terraflow has not been designed to output landscape decomposition into watersheds at a given scale/size but its result is used as an input for hierarchical watershed decomposition - see here http://terrain.cs.duke.edu/projects.php (and the papers by Arge et al. and Verdin and Verdin) However, the method used there can still leave you with small internal subwatersheds (r.watershed will have that too) that cannot be lumped with any bigger watershed. Helena
Working with Landsat Data
Background
- Landsat TM and ETM+ sensor information
- Data information and citation examples
- Pan-sharpening example in GRASS via Brovey transform
- Simple landscape classification example
Data Sources in order of ease of use
- University of Maryland's Global Land Cover Facility
- CaSiL Landsat Archive for CA
- USGS Global Viewer: free ortho-rectifed imagery
Sample GRASS commands
d.rgb r=landsat_band_3 g=landsat_band_2 b=landsat_band_1
#fuse bands into a single image:
r.composite r=landsat_band_3 g=landsat_band_2 b=landsat_band_1 out=landsat_bands_321
#blend landsat composite image with DEM data
#note that blending will require some tweaking: good results obtained with:
r.blend second=shaded_relief first=landsat_321 out=landsat_shade perc=6
#re-composite the blended results
r.composite landsat_shade.r g=landsat_shade.g b=landsat_shade.b out=landsat_shade.rgb
#create histogram equalized inputs (Fig 5)
r.colors map=landsat_band_1 color=grey.eq
r.colors map=landsat_band_2 color=grey.eq
r.colors map=landsat_band_3 color=grey.eq
r.composite r=landsat_band_3 g=landsat_band_2 b=landsat_band_1 out=landsat_eq_321
Simple comparison between the default output from r.composite on landsat bands 3,2,1 and the enhanced output from i.landsat.rgb.
 r.composite default
|
 i.landsat.rgb enhanced
|