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
