Hi guys
I've attached the code I have written for using OGR with Floatcanvas.
This lets you use geographic vector formats such as shapefiles.
I have not yet used the GDAL library (of which OGR is part of) for
raster formats. I must warn that the following code is in a very untidy
state with minimal comments as it is a work in progress, but feel free
to ask questions. Unfortunately it is part of a larger project, and so
will not run 'as is'. Hopefully you will be able to make it work (just
needs the wx window etc I think), otherwise I can probably put something
together.
To get OGR, apt-get install python-gdal, or go to their website
http://www.gdal.org/ogr/
ogr_shp.py contains a class with associated methods. This is used to
read the vector file and convert the data to a list of coordinates.
OGR supports spatial queries which is very handy (it also supports SQL
queries, but I havn't done anything with that).
The object projdat is for converting between various projections and
datums, so can probably be ignored at this stage.
MapCanvas.py is my class that contains the Floatcanvas (specifically
Navcanvas) and associated methods.
Theres not a whole lot thats important here - DrawOutlines and
DrawOutline are the ones we are interested in.
You might notice that I have done my own world to pixel coordinate
conversion, which is probably what led to my confusion with bounding
boxes. What is the most efficient and tidy way of doing this? I have
noticed conversion functions in Floatcanvas, I suppose this is what I
should be using? :)
Finally we need some code to tie the two together, which is fairly simple.
from ogr_shp import OGR_SHP
# Don't forget you will need to have a MapCanvas instance, or even just
a Floatcanvas will do with the right supporting functions
self.lat = -37.8
self.lon = 175.63
self.width = 10000
self.height = 10000
# Create an OGR_SHP object with the roads shapefile
roads = OGR_SHP('/home/sam/shapefiles/roads.shp')
# Apply a spatial filter using the lat, lon, and a buffer
# If you didn't want to apply a filter, you can just call
roads.get_line_features()
features = roads.spatial_filter_ll(self.lat,self.lon,self.width,self.height)
# nav_overview is an instance of MapCanvas
self.nav_overview.canvas.ClearAll()
self.nav_overview.features = {}
self.nav_overview.SetOrigin((self.lat,self.lon))
self.nav_overview.AddFeatures('roads',features, 'line')
self.nav_overview.DrawOutlines('roads')
All things going well, you should see the roads (or whatever) drawn on
the Floatcanvas. Of course, you will need a shapefile to run this with!!
'line' could be replaced with 'polygon' depending on the shapefile (ie
roads vs state boundaries).
Hope this was useful, any questions feel free to ask :)
Cheers, Sam
Christopher Barker wrote:
Sam,
Could you please join the list for further discussion -- I like to have
it there, and non-members can't post directly.
I'd be happy to send you some code, however its very messy at the
moment, so I'll send it once I've tidied it up a bit :)
OK, but I really am interested.
Turns out I passed in self explicitly by accident.
I thought so.
For anyone working with the bounding
box, it is important to be aware that you are working with pixel
coordinates.
I don't think so. I'm pretty sure the Bounding Box is in World coords.
This is quite purposeful -- all DrawObjects are defined in World Coords,
and pixel coords are only computed when drawing. So all the Bounding Box
checking, etc is done in World coords.
-Chris
# Thanks to Flavio Codeco Coelho for OGR code http://pyinsci.blogspot.com/2007_06_01_archive.html
import ogr, sys #, rpy
#import nzmg,nztm
from projdat import ProjDat
class OGR_SHP:
def __init__(self, filename):
self.map = None
query = None
self.projdat = ProjDat()
self.map = ogr.Open(filename)
print "Shapefile has %d layers" % self.map.GetLayerCount()
self.layer = self.map.GetLayer(0)
print "Layer 0 has %d features" % self.layer.GetFeatureCount()
def get_num_features(self):
return self.layer.GetFeatureCount()
def extract_double_field_from_feature(self, fieldname, feature):
feat = self.layer.GetFeature(feature)
index = feat.GetFieldIndex(fieldname)
return feat.GetFieldAsDouble(index)
def extract_string_field_from_feature(self, fieldname, feature):
feat = self.layer.GetFeature(feature)
index = feat.GetFieldIndex(fieldname)
return feat.GetFieldAsString(index)
def spatial_filter_bb(self,bb):
# Bounding box needs to be in geographic coordinates here
# So far, using NZMG
[[min_x,min_y],[max_x,max_y]] = bb
(min_y,min_x) = self.projdat.nztm2nzmg(min_y,min_x)
(max_y,max_x) = self.projdat.nztm2nzmg(max_y,max_x)
self.layer.SetSpatialFilterRect(min_x,min_y,max_x,max_y)
print "Spatial filter has %d features" % self.layer.GetFeatureCount()
#return self.get_features()
type = self.layer.GetLayerDefn().GetGeomType()
if type == 3: #polygons
return self.get_polygon_features()
if type == 2: #lines
return self.get_line_features()
def spatial_filter_ll(self, lat,lon,width,height):
# Using NZMG, should really check
(lly,llx) = self.projdat.geod2nzmg(lat,lon)
ury = lly + height
urx = llx + width
self.layer.SetSpatialFilterRect(llx,lly,urx,ury)
print "Spatial filter has %d features" % self.layer.GetFeatureCount()
#defn = self.layer.GetLayerDefn()
print 'layer type is', self.layer.GetLayerDefn().GetGeomType()
type = self.layer.GetLayerDefn().GetGeomType()
if type == 3:
return self.get_polygon_features()
if type == 2:
return self.get_line_features()
def get_line_features(self):
X = []
Y = []
features = []
feat = self.layer.GetNextFeature()
while feat:
geo = feat.GetGeometryRef()
for i in range(geo.GetPointCount()):
x = geo.GetX(i)
y = geo.GetY(i)
(y,x) = self.projdat.nzmg2nztm(y,x)
X.append(x)
Y.append(y)
feature = zip(X,Y)
X = []
Y = []
features.append(feature)
feat = self.layer.GetNextFeature()
#break
return features
def get_polygon_features(self):
X = []
Y = []
features = []
feat = self.layer.GetNextFeature()
while feat:
geo = feat.GetGeometryRef()
feature = []
if geo.GetGeometryCount()<2:
g1 = geo.GetGeometryRef( 0 )
for i in range(g1.GetPointCount()):
x = g1.GetX(i)
y = g1.GetY(i)
(y,x) = self.projdat.nzmg2nztm(y,x)
X.append(x)
Y.append(y)
feature = zip(X,Y)
X = []
Y = []
features.append(feature)
else:
for c in range( geo.GetGeometryCount()): # Whakatane has 2
ring = geo.GetGeometryRef ( c )
#print ring.GetPointCount()
for i in range(ring.GetPointCount()):
x = ring.GetX(i)
y = ring.GetY(i)
(y,x) = self.projdat.nzmg2nztm(y,x)
X.append(x)
Y.append(y)
feature = zip(X,Y)
X = []
Y = []
features.append(feature)
#X = []
#Y = []
#features.append(feature)
feat = self.layer.GetNextFeature()
#print feat
return features
import wx
import sys
import Image
from wx.lib.floatcanvas.Utilities import BBox
from wx.lib.floatcanvas import NavCanvas, Resources
from wx.lib.floatcanvas import FloatCanvas as FC
sys.path.append('../dem')
from projdat import ProjDat
from Node import NodeList
from shapely.geometry import *
class MapCanvas():
def __init__(self, panel,statusbar):
self.nav = NavCanvas.NavCanvas(panel, -1)
self.canvas = self.nav.Canvas
self.nodes = NodeList(self.canvas)
self.statusbar = statusbar
self.width = 0
self.height = 0
self.origin = (0,0)
self.world_bounding_box = None
self.pixel_bounding_box = None
self.projdat = ProjDat()
self.__bind()
def __bind(self):
# Bind mouse click for navcanvas
FC.EVT_LEFT_DOWN(self.nav, self.OnLeftDown )
# Bind mouse move for coords
#FC.EVT_MOTION(self.nav, self.OnMotion)
def AddBitmap(self, bitmap):
self.canvas.ClearAll()
buf = wx.Bitmap(bitmap,wx.BITMAP_TYPE_ANY)
(self.width,self.height) = buf.GetSize()
self.canvas.AddScaledBitmap(buf, (0,0), Height = self.height, Position = "bl" ) # top left
self.nav.ZoomToFit(self)
def AddFeatures(self,key,features,type):
self.features[key] = (type,features)
def DrawOutlines(self,key,colour='Black'): # rename function to something like DrawOutlinesFromFeatures
(type,features) = self.features[key]
for feature in features:
#print feature
self.DrawOutline(feature,type,colour)
if self.pixel_bounding_box == None:
self.nav.ZoomToFit(None)
else:
# Bounding box needs to be in pixels here
#print self.bounding_box
self.canvas.ZoomToBB(self.pixel_bounding_box)
def DrawOutline(self,feature,type,colour='Black'):
new_feature = []
(a,b) = self.nztm_origin
#print colour
for (x,y) in feature:
y = (y-a)/100
x = (x-b)/100
#nztm_y = float(a + y * 100)
#nztm_x = float(b + x * 100)
new_feature.append((x,y))
if type == 'polygon':
self.canvas.AddPolygon(new_feature, LineWidth = 1, LineColor = colour)#,FillColor = 'Red' , FillStyle = 'Solid')
if type == 'line':
self.canvas.AddLine(new_feature, LineWidth = 2, LineColor = colour)#,FillColor = 'Red' , FillStyle = 'Solid')
def SetOrigin(self,origin):
self.origin = origin
(lat,lon) = origin
self.nztm_origin = self.projdat.geod2nztm(float(lat), float(lon))
_______________________________________________
FloatCanvas mailing list
[email protected]
http://mail.mithis.com/cgi-bin/mailman/listinfo/floatcanvas
