Code

All code is Copyright © 2012, Gary Sherman

The following code samples used in The Geospatial Desktop are available:

clip_raster.sh
Shell script to clip a raster in GRASS using a vector layer

#!/bin/sh
############################################################################
#                          clip_raster.sh                                  #
#                          Copyright (C) 2012, Gary Sherman                #
# Clip a GRASS raster using a vector layer and mask                        #
#                                                                          #
############################################################################
 
# import the DRG
r.in.GDAL input=i61149e3_albers.tif output=ancc7_collars
# Extract the quad boundary from the boundary map
v.extract input=itma output=ancc7 where="TILE_NAME=’ANCC7’"
# Set the region to that extracted vector
g.region vect=ancc7
# Convert the extracted vector quad feature to a raster map
v.to.rast input=ancc7 output=ancc7_itma use=val
# Set the region to operate on to that of our DRG
g.region rast=ancc7_collars
# Set the mask for the operation to the raster created from the
# quad boundary vector
g.copy ancc7_itma,MASK
# Use map algebra to create the "clipped" raster
r.mapcalc ancc7=ancc7_collars
# Delete the mask
g.remove MASK
convert_shapefiles.sh
Convert all shapefiles in a directory to WGS 84 projection using ogr2ogr

#!/bin/bash
############################################################################
#                          convert_shapefiles.sh                           #
#                          Copyright (C) 2012, Gary Sherman                #
#           Convert all shapefiles in the current directory to             #
#           WGS84 projection. The converted shapefiles are placed          #
#           in the geo subdirectory.                                       #
#                                                                          #
############################################################################
 
for shp in *.shp
do
echo "Processing $shp"
ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:4326 geo/$shp $shp
done
create_hillshade.sh
Create a GRASS hillshade

#!/bin/sh
############################################################################
#                           create_hillshade.sh                            #
#                          Copyright (C) 2012, Gary Sherman                #
# Creates a GRASS hillshade raster from a DEM and shaded  relief raster    #
#                                                                          #
############################################################################
 
r.shaded.relief map=ancc6_new shadedmap=ancc6_shade zmult=4 --overwrite
cat myelevation.rules |r.colors map=ancc6_dem color=rules
r.his -n h_map=ancc6_new i_map=ancc6_shade r_map=ancc6_r \
   g_map=ancc6_g b_map=ancc6_b --overwrite
r.composite -d red=ancc6_r blue=ancc6_b green=ancc6_g \
   output=ancc6_combined --overwrite
import_volcanoes.py
Python script to parse a delimited text file of volcano data and create a shapefile

#! /usr/bin/env python
"""Parse a delimited text file of volcano data and create a shapefile.
   import_volcanoes.py                          
   Copyright (C) 2012, Gary Sherman                 
"""
 
# import the csv module
import csv 
# import the OGR modules
import osgeo.ogr as ogr 
import osgeo.osr as osr
 
# use a dictionary reader so we can access by field name
reader = csv.DictReader(open("volcano_data.txt","rb"), 
    delimiter='\t',     
    quoting=csv.QUOTE_NONE)
# END:imports
# START:ogr_setup
# set up the shapefile driver
driver = ogr.GetDriverByName("ESRI Shapefile") 
 
# create the data source
data_source = driver.CreateDataSource("volcanoes.shp") 
 
# create the spatial reference
srs = osr.SpatialReference() 
srs.ImportFromEPSG(4326)
 
# create the layer
layer = data_source.CreateLayer("volcanoes", srs, ogr.wkbPoint) 
 
# Add the fields we're interested in
field_name = ogr.FieldDefn("Name", ogr.OFTString) 
field_name.SetWidth(24) 
layer.CreateField(field_name) 
field_region = ogr.FieldDefn("Region", ogr.OFTString)
field_region.SetWidth(24)
layer.CreateField(field_region)
layer.CreateField(ogr.FieldDefn("Latitude", ogr.OFTReal)) 
layer.CreateField(ogr.FieldDefn("Longitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Elevation", ogr.OFTInteger)) 
# END:ogr_setup
# START:process_file
# Process the text file and add the attributes and features to the shapefile
for row in reader: 
  # create the feature
  feature = ogr.Feature(layer.GetLayerDefn()) 
  # Set the attributes using the values from the delimited text file
  feature.SetField("Name", row['Name']) 
  feature.SetField("Region", row['Region'])
  feature.SetField("Latitude", row['Latitude'])
  feature.SetField("Longitude", row['Longitude'])
  feature.SetField("Elevation", row['Elev']) 
 
  # create the WKT for the feature using Python string formatting
  wkt = "POINT(%f %f)" %  (float(row['Longitude']) , float(row['Latitude'])) 
 
  # Create the point from the Well Known Txt
  point = ogr.CreateGeometryFromWkt(wkt)  
 
  # Set the feature geometry using the point
  feature.SetGeometry(point) 
  # Create the feature in the layer (shapefile)
  layer.CreateFeature(feature) 
  # Destroy the feature to free resources
  feature.Destroy() 
 
# Destroy the data source to free resources
data_source.Destroy()
parse_earthquakes.rb
Ruby script to prep a text file of earthquake events with fixed length records to be imported as delimited text

#!/usr/local/bin/ruby
############################################################################
#                             parse_earthquakes.rb                         #
#                          Copyright (C) 2012, Gary Sherman                #
# Prep a text file of earthquake events with fixed length records to be    #
# imported as delimited text. The "|" is used as the delimiter.            #
#                                                                          #
############################################################################
 
f = File.open("../db_search10423")
# Skip the first two header records
2.times {f.gets}
# print the delimited header row containing the fields we are interested in
print "event_date|event_time|latitude|longitude|depth|magnitude\n"
# process the earthquake records
while not f.eof
  record = f.gets
  # use a fixed length approach to get the fields we want since
  # splitting on white space isn't feasible
  event_date = record[1..10]
  event_time = record[13..22]
  latitude = record[26..32]
  longitude = record[37..44]
  longitude_direction = record[46..46]
  depth = record[50..54]
  magnitude = record[66..69]
  # if the longitude is in the western hemisphere, it must be
  # negative
  if longitude_direction == 'W'
    longitude = -1 * longitude.to_f
end
  # print a delimited record
  STDOUT << event_date << "|" << event_time << "|" \
    << latitude << "|" << longitude << "|" \
    << depth.strip << "|" << magnitude.strip << "|\n"
end
# close the input file
f.close