Stuff collected by non-geo-scientists.
Widely used libraries:
Commonly used is WGS84.
Sometimes metric systems are useful. UTM zones are applicable worldwide, eg.
but the UTM system does not take into account country borders.
These CRS do
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
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))
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()
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
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