0

I have imported a point shapefile called site_locations.shp into a personal geodatabase. Now I'm trying to retrieve the data from the database using SQL:

dbpath = os.path.normpath(r'../geodata/site_locations.mdb') #<br>
msacon = pypyodbc.win_connect_mdb(dbpath) #<br>
msacur = msacon.cursor() #<br>
sqlstr = "SELECT ID, Shape, Notes FROM site_locations" #<br>
msacur.execute(sqlstr) #<br>
rows = msacur.fetchall() #<br>

I did the same with a sqlite database and I got the geometry in the so-called well known binary format (WKB). I was able to convert that WKB to well known text (WKT) and then to objects with x and y values. Apparently the field Shape does not contain data in WKB format, but in GeomSHP format. Is there a way in Python to convert this GeomSHP to WKB or so? I understand that some spatial databases have a function called CGeomWKB. Is there a Python implementation?

  • If you are working with ArcGIS 10.l you can use SearchCursor.da to retrieve a features shape in WKB. Just use the 'SHAPE@WKB' token. See here for more info. – JKyllo Apr 21 '15 at 20:55
1

Use GDAL/OGR to execute SQL on either the shapefile, or the ESRI Personal GeoDatabase (PGeo).

from osgeo import ogr
ogr.UseExceptions()
ds = ogr.Open(path_to_data_source)
lyr = ds.ExecuteSQL('SELECT ID, Notes FROM site_locations')
for idx in range(lyr.GetFeatureCount()):
    feat = lyr.GetFeature(idx)
    geom = feat.GetGeometryRef()
    wkb = geom.ExportToWkb()  # here it is
    feat = geom = None
ds = lyr = None  # close dataset

There are also other export formats, such as WKT, GML, KML and GeoJSON.

| improve this answer | |
  • Any solution in pure Python? – Dobedani Apr 22 '15 at 14:41
  • Nothing that I'm aware of, but you might be able to find some way to map the data structures from pyshp (which reads shapefiles) to geomet (which writes WKB). Both of these are pure Python libraries. – Mike T Apr 22 '15 at 19:36
1

Ok, I had a look at shapefile.py and was able to arrive at this solution:

from struct import unpack
points = []
msacur.execute("SELECT ID, Shape FROM site_locations")

# Field Shape contains binary data in GeomSHP format - to be converted
rows = msacur.fetchall()
for row in rows:
    # Find out the shape type
    shapeType = unpack("<i", row[1][0:4])[0] 
    if not shapeType in (1, 11, 21):
        print "Found shape is not a point - ignored ..."
    else:
        # Possible shape types are POINT, POINTZ and POINTM - 1, 11 and 21 resp.
        point = unpack("<2d", row[1][4:20])
        x = point[0]
        y = point[1]
        points.append({"id":int(row[0]), "lat":y, "lon":x})
        # TODO: parse remaining 8 or 16 bytes in case of POINTZ / POINTM 

Cheers!

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.