This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision | ||
spatialpython [2014/08/06 16:20] mantis [Between positions from a Cartesian coordinate system] |
spatialpython [2019/01/29 16:48] (current) mantis [Geocoding] |
||
---|---|---|---|
Line 1: | Line 1: | ||
+ | Stuff collected by non-geo-scientists. | ||
+ | |||
+ | Widely used libraries: | ||
+ | * apt-get install [[http://www.gdal.org/|python-gdal]] | ||
+ | * apt-get install python-geopy | ||
+ | * apt-get install python-pyproj # Cython wrapper to provide python interfaces to [[http://trac.osgeo.org/proj/|PROJ.4]] functions. | ||
+ | * apt-get install python-shapely | ||
+ | |||
+ | |||
+ | |||
+ | ======= Brush-up: spatial reference systems ======= | ||
+ | |||
+ | Commonly used is WGS84. | ||
+ | |||
+ | Sometimes **metric** systems are useful. [[http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system|UTM zones]] are applicable worldwide, eg. | ||
+ | * [[http://spatialreference.org/ref/epsg/23033/|23033]] for Austria, | ||
+ | * ... for France | ||
+ | * ... for Germany | ||
+ | * ... for Boston | ||
+ | |||
+ | |||
+ | but the UTM system does not take into account country borders. | ||
+ | |||
+ | These CRS do | ||
+ | |||
+ | * [[http://spatialreference.org/ref/epsg/31287/|Austria Lambert ]] for Austria | ||
+ | * ... for France | ||
+ | * ... for Germany | ||
+ | * ... for Massachusetts | ||
+ | |||
+ | |||
+ | ======= Convert coordinates from/to WGS84 ======= | ||
+ | |||
+ | <code bash> | ||
+ | sudo aptitude install python-pyproj | ||
+ | </code> | ||
+ | |||
+ | <code python> | ||
+ | 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 | ||
+ | |||
+ | </code> | ||
+ | |||
+ | ======= 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 [[http://code.google.com/p/geopy/downloads/list|download]] it here and install it in your python path). For a long introduction read the [[http://code.google.com/p/geopy/wiki/GettingStarted|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): | ||
+ | |||
+ | * Yahoo completely redesigned their services: http://developer.yahoo.com/boss/geo/ - 2000 queries per day are free. | ||
+ | * [[http://code.google.com/apis/maps/documentation/geocoding/| Google Maps API]] allows 2,500 geolocation requests per day. But at least no key needed here. | ||
+ | |||
+ | |||
+ | |||
+ | <code python> | ||
+ | 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)) | ||
+ | | ||
+ | | ||
+ | |||
+ | </code> | ||
+ | |||
+ | |||
+ | |||
+ | ======= ESRI shapefiles ======= | ||
+ | The Geospatial Data Abstraction Library (GDAL) is a C/C++ library capable of handling numerous raster and vector formats. Python bindings ([[http://www.pygis.de/index.php/GDAL|python-gdal]]) can be installed in Lucid: | ||
+ | |||
+ | <code bash> | ||
+ | sudo apt-get install python-gdal | ||
+ | </code> | ||
+ | Copying fields from an existing shapefile and adding new fields: | ||
+ | |||
+ | <code python> | ||
+ | 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() | ||
+ | </code> | ||
+ | |||
+ | |||
+ | ======= Aerial distance in meters ======= | ||
+ | |||
+ | ===== Between WGS84 positions ===== | ||
+ | |||
+ | |||
+ | From [[http://gis.stackexchange.com/questions/4022/looking-for-a-pythonic-way-to-calculate-the-length-of-a-wkt-linestring|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): | ||
+ | <code bash> | ||
+ | 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 | ||
+ | |||
+ | </code> | ||
+ | |||
+ | ===== Between positions from a Cartesian coordinate system ===== | ||
+ | |||
+ | |||
+ | <code python> | ||
+ | from shapely.geometry import Point | ||
+ | Point(0,0).distance(Point(1,1)) | ||
+ | </code> | ||
+ | |||
+ | |||
+ | Or if you already have wkt to work with: | ||
+ | |||
+ | <code python> | ||
+ | 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 | ||
+ | </code> | ||
+ | |||
+ | ======= Visualising maps ======= | ||
+ | |||
+ | ===== colourcode percentages ===== | ||
+ | |||