User Tools

Site Tools


spatialpython

Stuff collected by non-geo-scientists.

Widely used libraries:

  • apt-get install python-gdal
  • apt-get install python-geopy
  • apt-get install python-pyproj # Cython wrapper to provide python interfaces to PROJ.4 functions.
  • apt-get install python-shapely

Brush-up: spatial reference systems

Commonly used is WGS84.

Sometimes metric systems are useful. UTM zones are applicable worldwide, eg.

  • 23033 for Austria,
  • … for France
  • … for Germany
  • … for Boston

but the UTM system does not take into account country borders.

These CRS do

  • Austria Lambert for Austria
  • … for France
  • … for Germany
  • … for Massachusetts

Convert coordinates from/to WGS84

sudo aptitude install python-pyproj
from pyproj import Proj, transform
 
wgs84 = Proj(proj='latlong',datum='WGS84')
austria_lambert = Proj(init='epsg:31287')
 
 
x2, y2 = transform(wgs84, austria_lambert, x1, y1) # there
x3, y3 = transform(austria_lambert, wgs84,  x2, y2) # and back to WGS84 again

Geocoding

A frequent task, geocoding is done rather easily by whipping up a script and using the geopy library (apt-get install python-geopy, or download it here and install it in your python path). For a long introduction read the geopy wiki. Below is the fast-food variant:

Pick a geocoder (I tried only these two, but do check out the others and report back):

import sys
 
try:
    from geopy import geocoders
    from geopy.geocoders.google import GQueryError
except ImportError:
    print "You do not seem to have the python geopy module installed.\nCheck http://www.geopy.org/ for information about installation and usage."
 
 
geocoder = geocoders.Google( domain='maps.google.at')
#geocoder = geocoders.YahooPlaceFinder('consumer_key', 'consumer_secret')
 
address = "Wilhelmstr. 19, 1120 Wien, Austria"
 
try:   
    matched_address,(lat,lon) =  geocoder.geocode(address,exactly_one=False)[0] # use first result from result list
    print lat,lon
except GQueryError, qe:
    print("No match found for %s" % (address))

ESRI shapefiles

The Geospatial Data Abstraction Library (GDAL) is a C/C++ library capable of handling numerous raster and vector formats. Python bindings (python-gdal) can be installed in Lucid:

sudo apt-get install python-gdal

Copying fields from an existing shapefile and adding new fields:

import ogr
ogrdriver = ogr.GetDriverByName("ESRI Shapefile")
viennashp = ogrdriver.Open("./map/aut9__________nw.shp",0) # 0: readonly, 1: read-write
layer =  viennashp.GetLayerByName("aut9__________nw")
viennaout = viennaoutshp = ogrdriver.CreateDataSource(viennaout)
layerout = viennaoutshp.CreateLayer("nw",geom_type=ogr.wkbMultiLineString)
 
#create new fields
newField = ogr.FieldDefn()
newField.SetName("coverage")
newField.SetType(ogr.OFTInteger)
newField.SetWidth(5)
layerout.CreateField(newField)
 
#create roadId field by using the existing field definition from the source shapefile
fieldDef = layer.GetFeature(0).GetFieldDefnRef('ID')
layerout.CreateField(fieldDef)
 
outfeatureDef = layerout.GetLayerDefn()
road = layer.GetNextFeature() # if we are interested in Vienna's road graph
while road:
        outfeature = ogr.Feature(outfeatureDef)
        outfeature.SetGeometry(road.GetGeometryRef())
        outfeature.SetField("ID", road.GetFieldAsString("ID"))
        coverage = 1 # or whatever your magic wants it to be
        outfeature.SetField("coverage",coverage)
        layerout.CreateFeature(outfeature) # now create new feature
        road.Destroy()
        outfeature.Destroy()
        road = layer.GetNextFeature()
 
viennashp.Destroy() # clean up
viennaoutshp.Destroy()

Aerial distance in meters

Between WGS84 positions

From SO, but with a modification because with lat/lon it is always safest to use named parameters (otherwise lat and lon might get reversed and the earth implodes):

from geopy import distance,Point
 
 
# located in Northern Austria
lat0,lon0=48.2242843328, 16.1884416921
lat1,lon1=48.2377701258,  16.1900290977
 
# a number of other elipsoids are supported too
distance.VincentyDistance.ELLIPSOID = 'WGS-84'
 
a = Point(longitude=lon0,latitude=lat0)
b = Point(longitude=lon1,latitude=lat1)
 
print distance.distance(a,b).meters

Between positions from a Cartesian coordinate system

from shapely.geometry import Point
Point(0,0).distance(Point(1,1))

Or if you already have wkt to work with:

from shapely.wkt import loads
 
l=loads('LINESTRING(3.0 4.0, 3.1 4.1)') # creates a shapely.geometry.linestring.LineString
print l.length

Visualising maps

colourcode percentages

spatialpython.txt · Last modified: 2019/01/29 16:48 by mantis