2

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.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
George
  • 521
  • 4
  • 15
  • 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: http://gis.stackexchange.com/questions/193210/how-to-parse-through-a-esri-gdb-file-folder-after-opened-with-fiona 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

1 Answers1

3

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)])}
 ....
Community
  • 1
gene
  • 54,868
  • 3
  • 110
  • 187
  • 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
  • look at New in my answer – gene May 12 '16 at 19:04
  • Gene, I seem to still be crashing when using the next function, but I do get similar outputs for the prior functions in your new GDB example. Just so I have an understanding, does calling schema and crs step you through the GDB structure to put you in a position to start using next to call the next feature? – George May 12 '16 at 20:27
  • schema, crs or driver are only shapefile or GDB metadata, nothing to do with the records of the shapefile (next feature, next , ...) – gene May 12 '16 at 20:38
  • Ah, so your use of them (schema, crs, etc.) was merely to illustrate where in the GDB I am looking? – George May 12 '16 at 21:06
  • So digging around some more, it appears the next() command is causing Fiona's "collection.py" file (line 318) to crash... this is where it says "return next(self.iterator)" Did I set something up bad? Or is something wrong with my GDB file? – George May 12 '16 at 21:36
  • Gene, so I think I figured out what my problem is... I think the GDB file I'm using got corrupted migrating it from Windows to Linux. I'll have to request someone to physically copy it since the archive manager I used may be suspect in causing the problem I have. Although the next command would crash the program, the for loop just simply ran and exited with no errors. I believe I may have an empty GDB with metatags. – George May 12 '16 at 23:08
  • Gene, Just for closure's sake I got around my script crashing whenever doing a layer.next() call by enumerating a fiona.listlayers() call instead. For whatever reason, next() is hostile towards the ESRI-generated GDBs I've used. – George May 25 '16 at 16:59