Tag Archives: Open Source GIS

MongoDB and Shapely

18 Nov

I wanted to connect my mongoDB to Shapely so that I could create buffers and do some other things, like spit out my points and buffers as WKT. Here is what I came up with.

from pymongo import Connection,GEO2D
from shapely.geometry import asShape,Point


for q in db.b.find({“loc”:{“$near”:[-106,35]}}).limit(1):
point=Point(q[“loc”][0],q[“loc”][1])  #this is passing mongoDB to Shapely

point.x    #just checking I didn’t make a mistake


t.area  #checking. Result is 12.a lot of number



Now I can use mongoDB and find all the points in a bounding box or near another point or near a point where type=Coffee Shop, then send that to Shapely for a buffer or to measure distance or whatever, then spit it out as WKT. Saw a tutorial on using shapely with GDAL/OGR and Fiona. Goodle that and see how to send it all back to a shapefile.


Shapefile Creator Website Using Cherrypy

8 Nov

I showed how to create a shapefile using python, now I want to make the tool public with a web based front-end. My python web framework of choice is Cherrypy. -a minimal web framework.

My website will allow you to enter Longitude, Latitude and a comment. I am working on a text area where you can paste data. Here ie my website:

And when you submit, you can download a ZIP of your shapefile.

Pretty simple. I don’t really know what the use case is – if you have GIS and Long/Lat you would just import it directly. And if you don’t have GIS, why do you need a shapefile? But I guess this shows how you can perform GIS data functions through the web using python. Maybe connect some toolboxes and use Arcpy to allow people to run your models with their data – without giving your model and logic away.

My code is not pretty and no built for production – it is built to work. Here it is.

import cherrypy
import shapefile
from cherrypy.lib.static import serve_file
import zipfile

class SHP(object):

def index(self):
<BODY><FORM NAME=’shp’ action=’paul’>Latitude: <input type=’text’ name=’lat’><BR>Longitude: <input type=’text’ name=’long’>
<BR>Comment: <input type=’text’ name=’comment’><BR><input type=’submit’ Value=’SUBMIT’>

index.exposed = True

def paul(self,lat,long,comment):
w = shapefile.Writer(shapefile.POINT)
w.point(Flong, Flat)
prj = open(“Web.prj”, “w”)
epsg = ‘GEOGCS[“WGS 84”,’
epsg += ‘DATUM[“WGS_1984”,’
epsg += ‘SPHEROID[“WGS 84”,6378137,298.257223563]]’
epsg += ‘,PRIMEM[“Greenwich”,0],’
epsg += ‘UNIT[“degree”,0.0174532925199433]]’
zipped = zipfile.ZipFile(‘YourShapefile.zip’, mode=’w’)
return “Download Shapefile: <a href=’/get?filepath=C:\Documents and Settings\user\Desktop\YourShapefile.zip’>YourShapefile</a>”

paul.exposed = True

def get(self, filepath):
return serve_file(filepath, “application/x-download”, “attachment”)
get.exposed = True

cherrypy.config.update({‘server.socket_host’: ‘’,
‘server.socket_port’: 8000,


Shapefiles in Revit: Method II

25 Oct

Need to get shapefiles from your GIS in to Revit? I will tell you how using open source software. But first, Revit does not like coverages of over 20 miles so you may need to crop your shapefile.

Download Qgis – a free and open source GIS.

Open your shapefile in Qgis.

It will appear in the table of contents on the left. Right click on it and select ‘Save as..’ Choose AutoCAD DXF.

You now have a DXF of your shapefile.

Want to bring AutoCAD – or Revit – to GIS? Do everything in Reverse. But first activate your plugins.

Go to ‘Plugins’ and ‘Manage Plugins’. Check DXF2SHP.

Now click DXF2SHP and select your DXF. Choose a name and location to save your shapefile. Select the type. Probably polyline.

Here is a DXF loaded

Now your CAD drawing can be used in a GIS.

Python, Point Shapefile and Qgis: Part II

24 Oct

In my first post, I showed how to create a shapefile from an excel sheet – but it was hard coded by cell. Not very useful.  Here is the code to loop through the sheet and create shapefile, then send a text telling me the job finished. The last thing to do is render the shapefile – Mapnik would work for this.
Here is the code:

import shapefile
import xlrd
import smtplib
book = xlrd.open_workbook(“NAME OF XLS FILE.XLS”)
sh = book.sheet_by_index(0)
w = shapefile.Writer(shapefile.POINT)

def readRow():

for rownum in range(sh.nrows):
w.point(sh.cell_value(rowx=rownum, colx=0), sh.cell_value(rowx=rownum, colx=1))
w.record(sh.cell_value(rowx=rownum, colx=2),sh.cell_value(rowx=rownum, colx=3))
print ‘Reading Record: ‘+ str(rownum)

print ‘Complete…’


message = ‘Job Complete’
server = smtplib.SMTP( ‘smtp.gmail.com’, 587 )
server.login( ‘YOUREMAIL@gmail.com’, YOUR PASSWD’ )
server.sendmail( ‘YOUR PHONE #@mms.att.net’, ‘TO PHONE #@mms.att.net’, message )

The code opens the excel file formatted as: Lat, Long, Type, Address. Creates a point shapfile with fields for type and address.  Create a function that loops through the excel sheet. The Long, Lat, Type and Address are always in the same column so we only need to replace the row number of our cell and we do that by using the rownum variable in the loop.

While running, it will print out the current rownumber then print Complete.. when done. Lastly, It will send a text message saying that it has finished.

Simple.  Less than 25 lines of code. Set this to run as a task in Windows or as a Cron Job.

Python, Point Shapefile and Qgis

23 Oct

I have been getting data in a CSV and bringing it to Qgis using ‘Import Delimited Text.’  I could probably write a Python script in Qgis to do this automatically, but I only recently started reading the Qgis Python Cookbook. I have some experience with a Python Library for shapefiles already, so I will show how to create a point shapefile using this library.

Download Shapefile.py. All my maps are in WGS84, making this very easy. Open Notepad or your IDE and write the following:

import shapefile
w = shapefile.Writer(shapefile.POINT)
w.point(-106.505325, 35.069018)
w.point(-106.509378, 35.067741)
w.point(-106.527414, 35.108816)
w.record(‘res burg’,’central’)
w.record(‘auto burg’,’menaul’)
w.record(‘auto theft’,’eubank’)

I have three points in Albuquerque, with two fields: type and street. Running this file will create a shapfile without a PRJ. That’s alright, open the shapefile in Qgis. You should be asked for the projection. Choose WGS 84. It will display correctly. Now the points are on the map.

Right click on the shapefile in the table of contents and save as… Make sure the CRS is set to WGS 84.

Now when you open the shapefile you will not need to specify a projection.

Simple, right? Well, I hate hard coding the coordinates. I want to read from a file and everyone seems to send me data in Excel files so I downloaded XLRD and XLWT to read and write Excel files.

Here is how to add two points from an Excel Spreadsheet – you will want to run this as a loop to load all the data instead of hard coding it again. I will leave that up to you. (UPDATE: I can’t leave you hanging. Here is a new post with the code for looping). Here is the code to read the Excel file and write the shapefile.

book = xlrd.open_workbook(“test.xls”)
sh = book.sheet_by_index(0)
w = shapefile.Writer(shapefile.POINT)
w.point(sh.cell_value(rowx=0, colx=0), sh.cell_value(rowx=0, colx=1))
w.point(sh.cell_value(rowx=1, colx=0), sh.cell_value(rowx=1, colx=1))
w.record(sh.cell_value(rowx=0, colx=2),sh.cell_value(rowx=0, colx=3))
w.record(sh.cell_value(rowx=1, colx=2),sh.cell_value(rowx=1, colx=3))

Is this any easier than importing a CSV and creating the shapefile? No. But what it lets me do is schedule the Python file to run at intervals and update the shapefile even if I am at home sleeping. When graveyard shift comes in, they can see the latest data without me having to create it. It also allows me to add other functionality and python scripts.

I can filter the CSV by type or by location for example. Better yet, when a point in the shapefile exceeds some value, I can send a text to my phone – like this:

import shapefile
import xlrd
import smtplib
book = xlrd.open_workbook(“realtime.xls”)
sh = book.sheet_by_index(0)
w.record(sh.cell_value(rowx=1, colx=2))

x = sh.cell_value(rowx=1, colx=2)
if x > 100:
message = “At Point: “+ str(sh.cell_value(rowx=1, colx=0)) +”,” + str(sh.cell_value(rowx=1, colx=1)) +”Levels equal: “+str(sh.cell_value(rowx=1, colx=2))
server = smtplib.SMTP( “smtp.gmail.com”, 587 )
server.login( ‘YourAccount@gmail.com’, ‘Your Password’ )
server.sendmail( ‘From Phone Number@mms.att.net’, ‘ToPhoneNumber@mms.att.net’, message )

print ‘Levels are normal’

All the concatenation of Message =”….” will create a string that says: “At Point (-106.4567, 35.3445) Levels Equal: 125” but will fill in the numbers using the actual coordinates and value.

There are programs – like CrimeView – that do this sort of thing automatically. But the way I showed was free and can be modified as you wish. Trying modifying the code of a proprietary software solution….Good Luck and Happy Hacking.

Leaflet.js Interactivity II: Allow Users to Add Points

29 Jun

In Part I, I showed how to use buttons to add and remove layers, and how to capture popup text and display it in a <DIV> element. In this post I will show how to allow users to create their own points with popups. Here is an image of the page with no input yet.

First, I created a blank map with just a Cloudemade Tile.

Second, I created a <DIV> to hold my form:

LAT:<input type=”text” id=”lat”>
LONG:<input type=”text” id=”long”>
POPUP:<input type=”text” id=”popuptext” value=”POPUPTEXT”>
<button onclick=”createPoint()”>Create Point</button>

Lastly, I create the function in the form:

function createPoint()
var lat = document.getElementById(‘lat’).value;
var long = document.getElementById(‘long’).value;
var pop = document.getElementById(‘popuptext’).value;
var latint = parseInt(lat);
var longint = parseInt(long);
var mora = new L.LatLng(latint, longint),
MSD = new L.Marker(mora);

To see it working go to EducationalFacilityPlanning.com/UserAddPoints.html

Try typing in 35, -106 for Lat and Long

This needs more work, like a function to allow the user to click on the map and get a long,lat pair:

map.on(‘click’, onMapClick);
var popup = new L.Popup();

function onMapClick(e) {
var latlngStr = ‘(‘ + e.latlng.lat.toFixed(3) + ‘, ‘ + e.latlng.lng.toFixed(3) + ‘)’;

popup.setContent(“You clicked the map at ” + latlngStr);


But, more importantly, a function to save the points so the user can recreate their map.


Happy Hacking


Leaflet.js Interactivity

19 Jun

I have been playing with Leaflet and while I am not a JavaScript pro, I wanted to make my map more of an application. Two things I have started working on: Adding buttons and capturing popups.

In Leaflet, to add or remove a layer you only need to call map.addLayer(layerName) and map.removeLayer(layerName). I created a function


and then I created a hyperlink

<a href=”#” onClick=”addPoints()”>

I grabbed a button image and used it as the link. When you click the button, the map will add a layer.

Click Button to Add Layer

Red School Points are Added When Clicked.

The remove button does the same thing but uses the map.removeLayer() function.

I also wanted to capture popup info and display it on the webpage so I can display more data and do things with it. To display a popup in the webpage, I modified the popup function from Bryan McBride. Instead of:



I took the text in the setContent(), deleted the openPopup(), and used:

var container = document.getElementById(“display”);
container.innerHTML = “the text from setContent()”;

I have a DIV at the top of the page with an id=”display”. This is where the info from the popup will show up.

Top Line Will Open With PopUp Info

PopUp Displayed

Here is my full HTML file. It connects to my GeoServer so it will not work properly but you can see the code.