tommy's scratchpad

My previous post declared my plugging my nose and jumping into Python GIS stuff. Advice and info from Sean have been valuable in getting familiar with things, though I’ll be the first to admit I’m still green (writing Python like a Perl / C hack 🙂 ).

So here’s one of my first attempts at solving a real world problem: reprojecting a bunch of points from a CSV file:

#!/usr/bin/python

import sys
import mapscript

if (len(sys.argv) == 1):

	print sys.argv[0] + " <csvfile>"
	sys.exit(1)
else:
	projInObj  = mapscript.projectionObj("init=epsg:32619")
	projOutObj = mapscript.projectionObj("init=epsg:4326")

	print "id,geom,zone"
	f = open(sys.argv[1], 'r')
	for line in f:
		s = line.strip()
		k = s.split(",")
		wkt = "POINT(" + k[1] + " " + k[2] + ")"
		shape = mapscript.shapeObj.fromWKT(wkt) # man, I love WKT!!
		shape.project(projInObj, projOutObj)
		print k[0] + "," + shape.toWKT() + "," + k[3]
	f.close()

And that’s it! So now you can take the CSV output with the geometry encoded as WKT and hook it up to MapServer with some simple OGR OVF syntax in your mapfile:

CONNECTIONTYPE OGR
CONNECTION "<OGRVRTDataSource>
<OGRVRTLayer name='nb'>
<SrcDataSource>./nb.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<GeometryField encoding='WKT' field='geom'/>
</OGRVRTLayer>
</OGRVRTDataSource>"

Mind you, you might question why to reproject when MapServer can do this for you, but I digress. There’s probably more eloquent ways to do this than what’s been done above via mapscript, eh?