Differences

This shows you the differences between two versions of the page.

--- spatialpython [2014/08/06 16:20]
mantis [Between positions from a Cartesian coordinate system]
+++ spatialpython [2019/01/29 16:48]
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 =====

Small heaps of code

User Tools

Site Tools

Differences

Page Tools