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.

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

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

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('

Enter ISO2 code of city: ') print '

Select category of place to search in and around the city

' for index,place in show_first.items(): print '||%s|||||||%s'%(index,place) place_choice = int(raw_input('

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

Enter 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 '

---------------------%s--------------------

'%place for k,v in result.items(): print k,'

','==========================' for item in v: print item print '



' 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('

Enter ISO2 code of city: ') print '

Select category of place to search in and around the city

' for index,place in show_first.items(): print '||%s|||||||%s'%(index,place) place_choice = int(raw_input('

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

Enter 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 '

---------------------%s--------------------

'%place for k,v in result.items(): print k,'

','==========================' for item in v: print item print '



' 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.