Small heaps of code

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
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.Yahoo( "INSRRT APPLICATION ID")

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

Page Tools
