Geo @ ObjectGraph

Earth and Environmental Science Blog

Hey there! Thanks for dropping by my blog! Take a look around
and grab the RSS feed to stay updated. See you around!

Python and Geometry

Introduction

GDAL/OGR is a handy and heavy duty geometric library in Python. If you use any GIS systems, you often see “Python Console” which allows you to calculate complicated spatial computation. In GIS software, you can do this using your mouse (e.g. ArcGIS Toolbox and Spatial Analysis Extension); however, you need this technique if you need to customize or automate the process. I start using this library since I’m taking GIS course at Long Island University. In this blog entry, I would like to have an introductory tutorial how to handle geometry in Python using GDAL/OGR. If you do not have a Python setup, the easiest way is to open up Python Console window in GIS system. ArcGIS has python console but if you do not have access to ArcGIS, I recommend to use QGIS which is an open source software for free and it’s available for various platform. Since I’m using Mac, I’ll show those examples on QGIS below.

Download QGIS

Let’s Get Started

In QGIS (or your console) open console.

Plugins -> Python Console

Open Console Window

You see the prompt

>>

This is where you will type the script. The first line that you will type is:

import ogr

Then hit enter. Now, you are ready to use them.

Defining Geometric Objects

There are a couple of ways to represent geometric objects in the script. The first method is to build the object from Well-known Text (WKT) . I’m showing the example below:

line = ogr.CreateGeometryFromWkt('LINESTRING(5 6,20 30)')

This will create a “line” object which holds two points (5,6) and (20,30) on x-y axis.

The other way is to create an empty line object then add two points later.

line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint_2D(5,6)
line.AddPoint_2D(20,30)

You need 3 lines in this case. I’ll stick with the first way to build objects because it’s more intuitive.

Inspect Objects

You can simply use print statement and it will dump the components of the object.

print line

The output will be

LINESTRING (5 6,20 30)

In case if you mistype WKT notation, you will get None object. This indicates that the object creation was failed.

Let’s say, I forgot to add the last y-point (30) in order to create line2 object.

line2 = ogr.CreateGeometryFromWkt('LINESTRING(5 6,20 )')

When you say

print line2

The output returns

None

Creating Buffer

With buffer method, another object could be generated. Let’s create a point object first.

point = ogr.CreateGeometryFromWkt('POINT(15 20)')

Let’s generate circle. The first argument would be the distance of buffer. The circle is nothing but a bunch of points so you have to specify the number of points. Here, we use 30 points.

buff = point.Buffer(10,30)

This will create a buffer which has the distance 10.

Calculating Intersection between the line and buffer.

trimmed = buff.Intersection(line)

You will see the new object, “trimmed“.

print trimmed
LINESTRING (8.83253564459444 12.132057031351106,19.369569454978425 28.991311127965478)

This is trimmed line in the circle. The next section is an exercise to perform the same method and you will see completed figure at the end.

Exercise

In this exercise, you will use latitude and longitude as points. Create following:

  • line between following two locations
    • -95.98/45.03
    • -113.73/40.11
  • A buffer which has distance 5 and it consists of 30 points around the following location
    • -100.46/39.16

After creating objects, calculate the intersection between the line and buffer. Display the buffer information.

Extra credits:

You can export object as KML fragment. KML is just another xml which you can read in QGIS.

You can export to KML as below

point.ExportToKML()

Copy the output within the single quote and copy and paste inside following KML template (remember, this output is just fragment of KML so we need to complete it).

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Placemark>
<name>MyGeometry</name>
<!-- Replace this line with fragment of KML which is exported from python -->
</Placemark>
</kml>

Then save it as point.kml.

Do the same procedure for following objects. The file name should be: line.kml, buff.kml, trimmed.kml.

line.ExportToKML()
buff.ExportToKML()
trimmed.ExportToKML()

Go back to QGIS and add 4 KML files by clicking add button .

Right click on one of layers select “Zoom to Layer Extent”.

If you see following picture (or similar shape), it’s sucess!

Pink line represents the intersection. I increased the line size.

Here is the answer:

import ogr
line = ogr.CreateGeometryFromWkt('LINESTRING(-95.98 45.03,-113.73 40.11)')
point = ogr.CreateGeometryFromWkt('POINT(-100.46 39.16)')
buff = point.Buffer(5,30)
trimmed = line.Intersection(buff)
point.ExportToKML()
line.ExportToKML()
buff.ExportToKML()
trimmed.ExportToKML()

Download KML, Python Script and QGIS project file: Intersection Files (Zip File, 45KB)

Additional Reading and Projections

In order to apply projection for results, you need to use AssignSpatialReference method in SpatualReference class.

See following blog entry:

http://www.perrygeo.net/wordpress/?p=4

Closing Comment

I met a founder of NumPy library at New York Python meetup. NumPy is the core of GDAL/OGR library in Python and it handles all scientific calculations. It’s not only handling geometric calculations but it also handles other scientific calculations such as, FFT, Matrix, and so on. He has been developing the library since 90′s and he is also running a business to visualize and analyze scientific data using the library. If you interested in scientific calculation in Python, I definitely recommend to look at this library.

NumPy Website

Reference

 

You can leave a response, or trackback from your own site.

One Response to “Python and Geometry”

  1. gene says:

    You can use Shapely, better to do that
    http://pypi.python.org/pypi/Shapely/1.2.7

Leave a Reply