Geospatial Development with Python

1. Geospatial Development with Python

2. Writing a Program in Python • A program is a collection of Python statements, functions, etc. • A script is a list of commands which are run by a program written in a specific language such as Python • Scripts can be edited by text editors or preferably by a specialized editor that also allows the script to run • We can use Python command line in Windows • However, the use of the command line is limited (it does not support testing scripts) • So, instead we use the Python script editor (an IDE) that has a better interface than the command line mode

3. IDLE • We will also continue using the IDLE environment for Python scripting • IDLE is available with Python download • Using the interactive Python interpreter is great for writing small code • However, it is not meant to be saved! • For larger, more complex multi-line code, it is better to write it as a script (rather than using the interactive Python interpreter, e.g., IDLE)

4. Writing Script > File > New Window • This gives you a new ‘Untitled’ scripting window • There is no Python prompt in this window • Try: print (“Hello Python World”) in the window • Because the script is not interactive, nothing would happen • It needs to be run as a program (script) to do anything • So, let’s save it > File >Save As Save it as “hello.py” in: C:\Python27\ArcGIS10.1\My Scripts

5. Write and Run the hello.py Script • Type: print “Hello Python World” in IDLE; save it in the directory in your Python 27 download directory • After saving the script, run it in IDLE > Run > Run Module • We get the following window

6. Install PythonWin • For the textbook exercises you need to download PythonWineditor • Download it from: • http://www.cgl.ucsf.edu/Outreach/pc204/pythonwin.html • Starting from ArcGIS 10.x you can also use Python in ArcGIS for Desktop application

7. Scripting • Scripting in ArcGIS allows you to write a program to run a sequence of tools • Python is a primary, text-based scripting language in ArcGIS • Scripts can be used as a tool, by adding them to your toolbox or as stand alone programs • Scripts can handle complex tasks, like: • Looping through and inserting between records in a database • Error handling • Allowing applications (e.g., MS Excel and R statistical package) to interoperate and communicate • Scripts can be programmed to automatically run at a specific time

8. Running an ArcGIS tool in Python • Insert the book’s DVD in your optical drive, and install it by double clicking the esri executable application. • Download data and exercises in ESRIPress

9. Geospatial Data • Information about an object located on the surface of Earth using coordinates • e.g., location of an antenna, shape of a river, outline of a stratigraphic unit • Geospatial data associate information with a specific location • e.g., salinity or temperature of water in a well • The data are defined by a series of coordinates (lat, long) in addition to other information about the attributes of the spatial object (e.g., temperature, type, elevation, …) • Data also comes with metadata, about: • datum, projection, coordinate system, units of data, date file was created, … • These geospatial data are stored in a geodatabase • e.g., data about parks, Chattahoochee River, and roads in Atlanta • We can query the geodatabase and find: • What is the area of the Piedmont Park • How many bridges are there across the Chattahoochee • How fat is GSU from the Airport

10. GIS Data Formats • Raster format • Used to store bitmapped images, e.g., aerial photograph • e.g., DEM • In DEM, each pixel represents height of a point on Earth surface • Vector format • Represent spatial data with points, lines, and polygon • such as ESRI’s shapefile (e.g., Atlanta.shp) orSimple Features, which is an openGIS standard for storing point, line, and polygon, and their attributes

11. Tools • Tools such as PostGIS or ArcCatalog allow us to develop geodatabase and find answers to these questions • We can also use Web services, for example, to find the coordinate of each point, e.g., using: http://tinygeocoder.com • And then use an application to find the distance between the two points, using Python

12. Python Libraries for Geospatial Development • There are two popular Python libraries for raster and vector data: • GDAL (Geospatial Data Abstraction Library) • Is for raster-based geospatial data; aailable for download at: http://gdal.org • Windows user use Frank Warmerdam’s “FWTools” open source GIS binaries for Windows (32 bit); available at: http://fwtools.maptools.org # includes gdal and ogr • OGR • Is for vector-based geospatial data and is available for download at: http://gdal.org/ogr • NOTE: These are now merged together and are downloaded together under the common name GDAL

13. GDAL Data Model • GDAL’s data model includes: • Dataset, which hold the raster data in a collection of raster bands • Raster Band: represents a band, channel, or layer within an image (e.g., RGB image has red, green, and blue components of the image) • Raster size: specifies the width of the image in pixels and overall height of the image in lines • Georeferencing transform converts from x,y coordinates into georeferenced coordinates (on the surface of Earth) • Affine transformation mathematical formula allowing operations such as X offset, Y offset, X scale, Y scale, horizontal shear, vertical shear • Coordinate system includes the projection and datum • Metadataprovide additional information about the dataset

14. GDAL example Python code • Use GDAL to calculate the average of the height values contained in a DEM: from osgeo import gdal, gdalconst dataset = gdal.Open(”DEM.dat”) band = dataset.GetRAsterBand(1) fmt= “<“ + (“h” * band.XSize) totHeight = 0 for y in range (band.Ysize): scanline = band.ReadRaster (0, y, 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

15. OGR • Datasource: represents the file (e.g., a country), has many: • Layers: sets of related data, e.g., terrain, contour lines, roads layers • Spatial reference: specifies datum and projection • Attributes: additional metadata about the feature • Geometry: specifies the shape and location of the feature. Are recursive, i.e., can have sub-geometries

16. OGR example Python code import osgeo.ogr from osgeo import ogr shapefile = ogr.Open(“TM_WORLD_BORDERS-0.3.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()

17. Projections • We need to know which projection our data uses • Projection: deals with mapping geospatial position onto a 2D Cartesian plane with special cartographic algorithms • Also need to know the datum assumed by the data • Datum: model of the Earth shape • The good thins is that Python does the conversion of projections and datums into each other with its pyproj library. Pyproj is available at: https://code.google.com/p/pyproj/downloads/list

18. Projections • A process that maps the 3D Earth into a 2D map • There are hundreds of projection algorithms, but none is perfect; all have distortions to area, shape, length, angle, etc. • We need to choose one whose distortion does not matterfor our purpose • For our geo-spatial analysis (plotting or doing calculations), we need to know what type of projection was used to create the data!

19. Coordinate System • Projected Coordinate System • 3D Earth is projected onto the 2D Cartesian plane, and then puts points on that plane • The coordinates refer to a point on a 2D map that represent 3D Earth • Need to define an origin and map units • Unprojected Coordinate System • These directly and accurately refer to a point on the 3D Earth surface, but make it difficult to calculate distance • Latitude and longitude are examples of unprojectedcoor. Sys.

20. UTM • Both types of coordinate systems require a reference, • For unprojectedlat/long origin is at the Greenwich observatory, and zero latitude represents equator • Universal Transverse Mercator (UTM) coordinate system divides the Earth into 60 zones; each with a different map projection to minimize error • The reference point (origin) is at the intersection of the equator and the central meridian for each zone • False northing and false easting are added to the distance in m from this reference to avoid negative numbers

21. Types of Projection • There are 3 main groups: • Cylindrical • 3D Earth (imagine it is a small, lit spherical lantern) is put inside a cylinder, and projected (2D) mathematically onto the cylinder • The cylinder is then unwrapped • Conic • The 3D Earth is projected onto a cone surface • E.g., Albers Equal-Area Projection • Azimuthal • The 3D Earth is projected directly onto a flat surface • The Earth is seen as if we look at it from space

22. Datum

23. pyproj Library • The pyproj library includes two classes: • Proj and Geod • Proj • Is a transformation class that converts lat, long to cartographic native map x,y coordinates (default is m), and vice versa • The following instantiates a Proj object: projection = pyproj.Proj (proj=‘tmerc’, ellipse=‘WGS84’) • The transform() function converts one projection into another

24. pyproj Example Python Code import pyproj UTM_X = 565718.523517 # UTM coordinates UTM_Y = 3980998.9244 #Use two Proj objects to define the UTM and lat/long projections srcProj = pyproj.Proj (proj = “utm”, zone = “17”, ellips=“clrk66”, units=“m”) dstProj = pyproj.Proj (proj = ‘longlat’, ellips=‘WGS84’, datum=‘WGS84’) #transform the coordinates into latitude and longitude values 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)

25. Geod • Performs great circle distance (between two points on Earth) and angle calculations • The fwd() function starts from a point, and given an azimuth and distance, finds the ending point and returns the back azimuth (angle from the end to the beginning point; ccw back azimuth is negative, e.g., azimuth=30, backazimuth is -150) • The inv() function takes the coordinates of two points and returns the forward and backward azimuth

26. Geod example Python Code • #given a lat, long position of a point, uses a geod object to calculate another point 10 km northeast of the first point angle = 15 #15 degrees from N distance = 10000 geod – pyproj.Geod(ellips=‘clrk66’) long2, lat2, invAngle = geod.fwd(long, lat, angle, distance) print “%o.4f, %0.4f is 10 km northeast of %0.4f, %0.4f” %(lat2, long2, lat, long)

27. Shapely Packagehttp://trac.gispython.org/lab/wiki/Shapely • Shapely is a python library for working and analyzing 2D geospatial geometries in the Cartesian plane 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)]) intesect = circle.intersection(square) for x,y in instersect.exterior.coords: print x, y

28. Example # location of two points given by their coordinates # It is best of the locations are stored in a geodatabase Import pyproj lat1, long1 = (37.37.775196, 122.419204809662) lat2, long2 = (33.790348, -117.226085) geod = pyproj.Geod (ellips =“WGS84”) # WGS84 is a Datum angle1, angle2, distance = geod.inv (long1, lat1, long2, lat2) Print “Distance is %0.2f m” % distance Distance is 2937.41 m

29. Mapping Applications • Mapnik takes geospatial data in any format (supported by GDAL/OGR) and renders them in a map • It provides rules and symbolizers to control how objects appear on the map • Mapnik supports Python

30. 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=“TM_WORLD_BORDERS-0.3.shp”) layer.style.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”) Mapnik is available from: http://mapnik.org

31. Work with geospatial data with Python • Before anything, download the following: • Make sure the download is for the version of the Python which is on your computer (e.g., v. 2.7) and for your machine (32 bit vs. 64 bit) • GDA/OGR http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries • pyprojhttps://code.google.com/p/pyproj/downloads/list • Shapely https://github.com/sgillies/shapely