0

After I successfully open up an ESRI GDB file/folder with Fiona, what open-source Python package should I use to do a spatial join?

I was told Shapely is not the appropriate package for this.

I am completely new to GIS.

  • 1
    More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc with Fiona and Shapely – gene May 12 '16 at 16:33
  • As long as you are specifying that you do not want to use Shapely I think this question is unique from the one linked by gene. However, why were you told Shapely was not appropriate for this application? – Evil Genius May 12 '16 at 17:28
  • I thought using Shapely was the most straightfoward way to approach doing a spatial join. I was told in this thread: gis.stackexchange.com/questions/193210/… that: "Shapely is a computational geometry library: buffers, intersections, &c. It doesn't do spatial joins but can be used to help your own join implementations." – George May 12 '16 at 18:41
2

Look first at More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc

The results of Fiona are Python dictionaries (GeoJSON format)

import fiona
layer = fiona.open("test_regex.shp")
# first feature
first = layer.next()
print first
{'geometry': {'type': 'Polygon', 'coordinates': [[(203371.23902535878, 89863.381050732), (203353.45178501407, 89474.60279748365), (203217.99038220354, 89246.7813473023), (203147.8656364849, 89284.24525254924), (203144.05848558623, 89691.61039870887), (203371.23902535878, 89863.381050732)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1), (u'dip', 30), (u'dip_dir', 150)])

Now if you want to use topological predicates (contents(), intersect() and others) the easier is to use Shapely (or equivalents, look at Python Geo_interface applications)

from shapely.geometry import polygon, shape
geom = shape(first['geometry'])
print geom
POLYGON ((203371.2390253588 89863.381050732, 203353.4517850141 89474.60279748365, 203217.9903822035 89246.7813473023, 203147.8656364849 89284.24525254924, 203144.0584855862 89691.61039870887, 203371.2390253588 89863.381050732))

Now you can use geom.intersection(another_geom), geom.contains(another_geom) and others therefore Shapely is an appropriate package for this (try to do the same with the original data)

You can also use Spatial Indexes (Fastest way to join many points to many polygons in python) and other solutions as GeoPandas, based on Fiona, Shapely and Pandas (More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2))

New

With a GDB file

gdb = fiona.open("test.gdb")
gdb.driver
u'FileGDB'
gdb.schema
{'geometry': 'MultiLineString', 'properties': OrderedDict([(u'id', 'int')])}
gdb.crs
{'init': u'epsg:31370'}
# first feature
gdb.next()
{'geometry': {'type': 'MultiLineString', 'coordinates':   [[(-1.2543001174926758, 0.22390007972717285), (-1.05430006980896, 0.6630001068115234), (-0.6935000419616699, 0.6284000873565674), (-0.30660009384155273, 0.7262001037597656), (0.3064999580383301, 0.8891000747680664)]]}, 'type': 'Feature', 'id': '1', 'properties': OrderedDict([(u'id', None)])}
# all the features
for feature in gdb:
   print feature
{'geometry': {'type': 'MultiLineString', 'coordinates': [[(-1.2543001174926758, 0.22390007972717285), (-1.05430006980896, 0.6630001068115234), (-0.6935000419616699, 0.6284000873565674), (-0.30660009384155273, 0.7262001037597656), (0.3064999580383301, 0.8891000747680664)]]}, 'type': 'Feature', 'id': '1', 'properties': OrderedDict([(u'id', None)])}
{'geometry': {'type': 'MultiLineString', 'coordinates': [[(-0.3847999572753906, 0.14569997787475586), (-0.16669988632202148, -0.12890005111694336), (0.09349989891052246, -0.3803999423980713)]]}, 'type': 'Feature', 'id': '2', 'properties': OrderedDict([(u'id', None)])}
 ....
  • Thank you. I just tried your first example, but my toy app crashes when I try the layer.next() command. I know I've opened up the GDB file "correctly" as Fiona does not provide an exception and I can see that there is metadata tags. – George May 12 '16 at 18:25
  • You have a problem because the result of fiona.open is a Python iterator with next() and if it crashes you have not opened up the GDB file "correctly" – gene May 12 '16 at 18:40
  • The Fiona open command doesn't cause an exception and I was told that I could do something like this the code below to verify successfully opening a GDB file with fiona.open(gdb_path) as src: print(src.meta) – George May 12 '16 at 18:45
  • Perhaps but if next() does not work you cannot access the features, so.. – gene May 12 '16 at 18:54
  • Is there any other way for me to establish / verify that I did open the GDB? I know I can print layer the first time it's assigned when Fiona opens up the GDB file. – George May 12 '16 at 18:59

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.