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 =====