Geospatial developement simplified with python

 

 

 

 

In this article I am going to show how easy it is to build a real world Geo-spatial application with python.

GIS?

What is Geospatial development or Geographical Information Systems (GIS)?. According to Wikipedia ,

Geographical information systems (GIS), which is a large domain that provides a variety of capabilities designed to capture, store, manipulate, analyze, manage, and present all types of geographical data, and utilizes geospatial analysis in a variety of contexts, operations and applications.

Basic Applications

Geo-spatial analysis, using GIS, was developed for problems in the environmental and life sciences, in particular ecology, geology and epidemiology. It has extended to almost all industries including defense, intelligence, utilities, Natural Resources (i.e. Oil and Gas, Forestry etc.), social sciences, medicine and Public Safety (i.e. emergency management and criminology), disaster risk reduction and management (DRRM), and climate change adaptation (CCA). Spatial statistics typically result primarily from observation rather than experimentation.

For more visit this page http://en.wikipedia.org/wiki/Geospatial_analysis

What Python got do with it?

If somebody asks you to find the locations of churches around 30KM from your new residence,as a first choice you go for web and find manually by observing the map of your city.But it will be complex ,if you wish to know about a location of many types. GIS are the systems those analyze earth data and fetches you the results.There are two types of earth data.

1)Raster data,data consisting scanned images of earth

2)Vector data,data consisting geometry drawn images

Analyzing both types of data is essential for many organizations like Military,Survey and few more.In initial days the people from mathematics background are hired for performing Geo-spatial development because all the basic operations,algorithms applied should be mathematically designed peer-peer.But now a days abstract bindings for existing libraries enables normal software developers to use wrappers for creating Geo-spatial applications.If you chose the correct programming language like Python,we can create powerful Geo-spatial applications in less time with less effort.Just programmer need to be familiar with terminology of geography like latitude,longitude,hemisphere,datum,meridian,directions,shape files etc.

Python got 3 fantastic binding libraries which are wrote upon existing libraries in c++

1) GDAL/OGR

2)PyQgis

3)  ArcPy

In this post i am going to show you a jump start example for finding the locations of a city in United States.I have United States city data with me.If your own region has data available,you can use the same code for finding theaters,parks etc around few Kilometers. My application is totally off-line,since I am analyzing a downloaded dataset.

I am going to use GDAL/OGR library for creating the application,You need to install GDAL/OGR, pyproj , shapely python libraries before starting to build application.For installing those libraries see

http://gis.stackexchange.com/questions/9553/whats-the-easiest-way-to-install-gdal-and-ogr-for-python/124751#124751

Seeing is believing

See the application running,and be confident.I am going to find cliffs around Texas,with 50 KM range

 

Here my main program is cliffs.py and program prompts to enter ISO2 code of city,TX-Texas CA-Callifornia.

 

Selection_017

now my application finally returns me all the cliffs around Texas in 50KM range

Selection_018

How I did that?

First we need to download required Datasets into our directory,two datasets we required here are

1) https://app.box.com/s/2s7i5culjrkq6sm3fqr5

2) https://app.box.com/s/7qzk2y3bpvo264hgxs0r

 

First file is Place-Find.zip,unzip it and save all files in the same directory of the program cliffs.py.Next second file is NationalFile_20141005.txt.Keep this file in same directory.Now we got the required datasets.

It will then look like

Selection_020

I am creating a new file called settings.py which stores the information of places to display

#settings.py
cats= []
with open('NationalFile_20141005.txt','r') as fil:
    for i in fil.readlines():
        cats.append((i.rstrip().split('|'))[2])

#List the places and take input from User for Park,Bar,Hotel etc
#we can use dict(enumerate(set(cats))) here but we need to delete FEATURE_CLASS,Unknown fields from list
show_first = {k:v for k,v in enumerate(set(cats)) if v != 'FEATURE_CLASS' and v != 'Unknown'} 

Once observe the  NationalFile_20141005.txt,it consists of information about the locations like parks,bays,cliffs etc.
Next I am going to create the main program cliff.py.Here go imports
from __future__ import division
from settings import show_first

from osgeo import ogr
import shapely.geometry
import shapely.wkt
osgeo deals with opening shape files,shapely is helpful in translating them into geometric shapes.division is used to convert KM into angular distance to measure(1 degree=100KM). We are importing show_first dictionary from settings.py
shapefile = ogr.Open("tl_2014_us_cbsa.shp")
layer = shapefile.GetLayer(0)

 this opens a file in the same directory and creates a shapefile object.Next we created a layer using GetLayer function,Since that shape file t1_2014-us_cbsa.shp has only one layer we can use GetLayer(0) to fetch that single layer.
city = raw_input('\nEnter ISO2 code of city: ')
print '\nSelect category of place to search in and around the city\n'
for index,place in show_first.items():
    print '||%s|||||||%s'%(index,place)
place_choice = int(raw_input('\nEnter code of place from above listing: '))
place = show_first[place_choice]
distance = int(raw_input('\nEnter with in range distance(KM) to find %s: '%place))
#converting distance range to angular distance $$$$ 100 KM = 1 Degree $$$
MAX_DISTANCE = distance/100 # Angular distance; approx 10 km.

Above lines are normal python code for intaking the preferences from user and last line is converting distance range to angular.

print "Loading urban areas..."

# Maps area name to Shapely polygon.
urbanAreas = {} 

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    shape = shapely.wkt.loads(geometry.ExportToWkt())
    dilatedShape = shape.buffer(MAX_DISTANCE)
    urbanAreas[name] = dilatedShape
In the above code we are translating geometry of each feature from the shape file into dilated shape and storing it in urbanAreas.
These urbanAreas now hold data for city polygons.Now i will open the NationalFile_20141005.txt and read the latitude and longitude values of the points around MAX_DISTANCE.If that point is within the polygon save it,else ignore.
f = open("NationalFile_20141005.txt", "r")
result = {}
for line in f.readlines():
    chunks = line.rstrip().split("|")
    if chunks[2] == place and chunks[3] == city:
        cliffName = chunks[1]
        latitude = float(chunks[9])
        longitude = float(chunks[10])
        pt = shapely.geometry.Point(longitude, latitude)
        for urbanName,urbanArea in urbanAreas.items():
            if urbanArea.contains(pt):
                if not result.has_key(cliffName):
                    result[cliffName]=[urbanName]
                else:
                    result[cliffName].append(urbanName)
 result dictionary is for saving the cliff_name,list of locations.The code above is quite obvious.Now print the results
print '\n---------------------%s--------------------\n'%place
for k,v in result.items():
    print k,'\n','=========================='
    for item in v:
        print item
    print '\n\n'
f.close()
this finishes our cliff.py.This application works entirely offline.Datasets are huge because of the details they are holding.The final code looks this way
#cliffs.py                    
from __future__ import division
from settings import show_first

from osgeo import ogr
import shapely.geometry
import shapely.wkt



shapefile = ogr.Open("tl_2014_us_cbsa.shp")
layer = shapefile.GetLayer(0)
#code for displaying and as

city = raw_input('\nEnter ISO2 code of city: ')

print '\nSelect category of place to search in and around the city\n'
for index,place in show_first.items():
    print '||%s|||||||%s'%(index,place)

place_choice = int(raw_input('\nEnter code of place from above listing: '))
place = show_first[place_choice]

distance = int(raw_input('\nEnter with in range distance(KM) to find %s: '%place))

#converting distance range to angular distance $$$$ 100 KM = 1 Degree $$$
MAX_DISTANCE = distance/100 # Angular distance; approx 10 km.
print "Loading urban areas..."

# Maps area name to Shapely polygon.
urbanAreas = {} 

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    shape = shapely.wkt.loads(geometry.ExportToWkt())
    dilatedShape = shape.buffer(MAX_DISTANCE)
    urbanAreas[name] = dilatedShape

print "Checking %ss..."%place


f = open("NationalFile_20141005.txt", "r")
result = {}
for line in f.readlines():
    chunks = line.rstrip().split("|")
    if chunks[2] == place and chunks[3] == city:
        parkName = chunks[1]
        latitude = float(chunks[9])
        longitude = float(chunks[10])
        pt = shapely.geometry.Point(longitude, latitude)
        for urbanName,urbanArea in urbanAreas.items():
            if urbanArea.contains(pt):
                if not result.has_key(parkName):
                    result[parkName]=[urbanName]
                else:
                    result[parkName].append(urbanName)


print '\n---------------------%s--------------------\n'%place
for k,v in result.items():
    print k,'\n','=========================='
    for item in v:
        print item
    print '\n\n'
f.close()
 This completes our application.We can do many other things like finding borders of countries accurately,length from one place to another,satellite functions.Any average python developer can excel in Geo-spatial informatics because he got powerful tool as programming language and open source libraries built already with hiding complexity behind great mathematical functions.code forthis application is available at GITHUB.
If you haven’t done programming in GIS,then it might look complex stuff to  you.But learn basics and see it once again,you feel this post nerdy.

6 thoughts on “Geospatial developement simplified with python

  1. If you’re looking to make a geospatial website, flask-admin just updated with geoalchemy2 support. It’s not very advanced yet, but it does allow you to edit individual PostGIS records.

  2. Great post. Good simple code, that cuts to the chase.
    It is worth mentioning that Neither arcpy nor pyqgis are used in this post, but both are amazing tools that far outstripe the capabilities of shapely. Shapely, however is incredible at what it does, and has a much lower learning curve.

    To use arcpy you need ArcGIS software installed on your machine, which has a costly license agreement and only works on windows (and they have a more expensive server version that also works on Linux).
    Pyqgis requires you to have QGIS software installed on your machine, it is free and open source and works on windows, mac osx, and a wide variety of linux distributions.
    In either case be prepared to spend a good bit of time reading the documentation.

  3. I’ll right away grab your rss feed as I can’t to find your email subscription hyperlink or newsletter service.
    Do you’ve any? Please permit me understand in order that I may
    just subscribe. Thanks.

  4. The GeoDjango tutorial shows a detailed GDAL/OGR usage. GDAL/OGR reads GIS “shape file” data into a data base (in this case PostGIS for PostgreSQL). I recommend that as good start to learn GIS.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s