那些年 Python 攻佔了 GIS / The Year Python Takes Over GIS

那些年Python攻佔GIS
The Year Python Takes Over GIS

鄧東波 Dongpo Deng
中央研究院資訊科學研究所
﹠
胡崇偉 marr, Tsung Wei Hu
中央研究院人文社會科學中心
PyCon 2012

Saturday, June 9, 2012

Content

• Introduction - GIS and why GIS uses
python
• Python-blinding core geospatial libraries
• the use of python in Desktop GISs
• Web application framework for geospatial


What is GIS?
• GIS is stand for Geographic
Information System
• integrates hardware, software,
and data for capturing,
managing, analyzing, and
displaying geospatial data.
• allows people to use methods
for understanding,
interpreting, and visualizing
relationships and patterns of
geospatial data

Why Geospatial domain uses Python

• Easy to learn
• Code is readable
• Large community
• Easy interaction with C and Java libraries
• Many existing modules and packages
• core geospatial libraries
• map rendering
• database
• web server

Picture from http://pypi.python.org/pypi/collective.geo.bundle


Geospatial development tasks
• Visualizing geospatial data


Geospatial development tasks
• Analyzing geospatial data
• e.g. How many people should escape as
Kuosheng nuclear power plant (核二廠)
explodes?
• Create a geospatial mashup


The geospatial development tasks
involve
• Math- analytic geometry
• e.g.Euclidean geometry
• Computer graphics (rendering)
• e.g. rending, visualizing
• Database
• General Search Tree (GiST)
• open geospatial standards,
• e.g. GML, WKT


Python libraries for geospatial
development
• Reading/ Writing geospatial data
• GDAL/OGR
• Dealing with Projections
• pyproj
• Analyzing and manipulating geospatial data
• Shapely
• Visualizing geospatial data
• Mapnik


GDAL
(Geospatial Data Abstraction Library)
• read through it one scanline at a time
from osgeo import gdal,gdalconst from ﬁle
import struct
dataset = gdal.Open("DEM.dat") • use the struct standard Python library
band = dataset.GetRasterBand(1) module to read the individual values out
fmt = "<" + ("h" * band.XSize) of the scanline.
totHeight = 0
for y in range(band.YSize): • Each value corresponds to the height of
scanline = band.ReadRaster(0, y, that point, in meters
band.XSize, 1, band.XSize,
1,band.DataType)
values = struct.unpack(fmt,
scanline)
for value in values:
totHeight = totHeight + value
average = totHeight / (band.XSize *
band.YSize)
print "Average height =", average

Source from: Westra, 2010, Python Geospatial Development

OGR
(OpenGIS Simple Features Reference Implementation)

• uses OGR to read through the contents of a Shapeﬁle,
• printing out the value of the NAME attribute for each
feature, along with the geometry type

from osgeo import ogr
shapefile = ogr.Open("TM_WORLD.shp")
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
name = feature.GetField("NAME")
geometry = feature.GetGeometryRef()
print i, name, geometry.GetGeometryName()


PyProj
• a location speciﬁed using UTM zone 17 coordinates.
• Using two Proj objects to deﬁne the UTM Zone 17 and
lat/long projections
• translates this location's coordinates into latitude and
longitude values
import pyproj
UTM_X = 565718.523517
UTM_Y = 3980998.9244
srcProj = pyproj.Proj(proj="utm", zone="11",
ellps="clrk66", units="m")
dstProj = pyproj.Proj(proj='longlat', ellps='WGS84',
datum='WGS84')
long,lat = pyproj.transform(srcProj, dstProj, UTM_X, UTM_Y)
print "UTM zone 17 coordinate (%0.4f, %0.4f) = %0.4f, %0.4f" %
(UTM_X, UTM_Y, lat, long)

Shapely GEOS
import shapely.geometry
pt = shapely.geometry.Point(0, 0)
circle = pt.buffer(1.0)
square = shapely.geometry.Polygon([(0, 0),
(1, 0),
(1, 1),
(0, 1),
(0, 0)])
intersect = circle.intersection(square)
for x,y in intersect.exterior.coords:
print x,y


Mapnik

• Mapnik is an open source mapping
toolkit for desktop- and server-based
map rendering, written in C++.
• there are Python bindings to facilitate
fast-paced agile development.
• OpenStreetMap project (OSM) uses
Mapnik in combination with an
Apache Web Server module
(mod_tile) to render tiles that make
up the OSM 'Slippy Map' Layer.



Mapnik
import mapnik
symbolizer = mapnik.PolygonSymbolizer(
mapnik.Color("darkgreen"))
rule = mapnik.Rule()
rule.symbols.append(symbolizer)
style = mapnik.Style()
style.rules.append(rule)
layer = mapnik.Layer("mapLayer")
layer.datasource = mapnik.Shapefile(
file="TW_village.shp")
layer.styles.append("mapStyle")
map = mapnik.Map(800, 400)
map.background = mapnik.Color("steelblue")
map.append_style("mapStyle", style)
map.layers.append(layer)
map.zoom_all()
mapnik.render_to_file(map, "map.png", "png")

Desktop GIS

Pic from http://www.pressebox.de/attachments/details/39739


Script Languages in ESRI family


ArcPy
• ArcPy is a package for performing geographic data
analysis, data conversion, data management, and
map automation in ArcGIS with Python.
• ArcPy includes three modules:
• mapping module (arcpy.mapping)
• Spatial analyst module (arcpy.sa)
• Geostatistical analyst module (arcpy.ga)

ArcGIS ArcPy Python

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/What_is_ArcPy/

arcpy in ArcGIS 10


Python in QGIS
• There’s a number of ways to access the layers in
QGIS.
• Each way starts by ﬁrst referencing the
QgisInterface class which is called iface in the
Python bindings.


Example


GeoDjango


http://140.109.160.129:8000/admin/world/worldborders/

Remarks

• There are many Python libraries or
applications for geospatial purposes
• Python is increasing its value in geospatial
domain
• Will Python take over GIS? .....Let’s see!


Thank you for attention!

dongpo.deng@gmail.com
marr.tw@gmail.com


那些年 Python 攻佔了 GIS / The Year Python Takes Over GIS

More Related Content

Similar to 那些年 Python 攻佔了 GIS / The Year Python Takes Over GIS

More from pycontw

Recently uploaded

那些年 Python 攻佔了 GIS / The Year Python Takes Over GIS