I have been trying to create a function that does basically the same thing that the QGIS "dissolve" function. I thought it would be super easy but well apparently not. So from what I have gathered around, the use of fiona with shapely should be the best option here. I just began to mess around with vector files so this world is pretty new to me and python as well.
For these example, I am working with a county shapefile founded here http://tinyurl.com/odfbanuSo here are some piece of code I gathered but can't find a way to make them work together
For now my best method is the following based on : http://sgillies.net/blog/870/a-more-perfect-union-continued/. It works fine and I get a list of the 52 states as Shapely geometry. Please feel free to comment if there is a more straight forward way to do this part.
from osgeo import ogrfrom shapely.wkb import loadsfrom numpy import asarrayfrom shapely.ops import cascaded_unionds = ogr.Open('counties.shp')layer = ds.GetLayer(0)#create a list of unique states identifier to be able#to loop through them laterSTATEFP_list = []for i in range(0 , layer.GetFeatureCount()) : feature = layer.GetFeature(i) statefp = feature.GetField('STATEFP') STATEFP_list.append(statefp)STATEFP_list = set(STATEFP_list)#Create a list of merged polygons = states #to be written to filepolygons = []#do the actual dissolving based on STATEFP#and append polygonsfor i in STATEFP_list : county_to_merge = [] layer.SetAttributeFilter("STATEFP = '%s'" %i ) #I am not too sure why "while 1" but it works while 1: f = layer.GetNextFeature() if f is None: break g = f.geometry() county_to_merge.append(loads(g.ExportToWkb())) u = cascaded_union(county_to_merge) polygons.append(u)#And now I am totally stuck, I have no idea how to write #this list of shapely geometry into a shapefile using the#same properties that my source.So the writing is really not straight forward from what I have seen, I really just want the same shapefile only with the country dissolve into states, I don't even need much of the attribute table but I am curious to see how you can pass it on from the source to the new created shapefile.
I found many pieces of code for writing with fiona but I am never able to make it work with my data.Example from How to write Shapely geometries to shapefiles? :
from shapely.geometry import mapping, Polygonimport fiona# Here's an example Shapely geometrypoly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])# Define a polygon feature geometry with one attributeschema = { 'geometry': 'Polygon', 'properties': {'id': 'int'},}# Write a new Shapefilewith fiona.open('my_shp2.shp', 'w', 'ESRI Shapefile', schema) as c: ## If there are multiple geometries, put the "for" loop here c.write({ 'geometry': mapping(poly), 'properties': {'id': 123}, })The problem here is how to do the same with a list of geometry and how to recreate the same properties than the source.Thank you!
أكثر...
For these example, I am working with a county shapefile founded here http://tinyurl.com/odfbanuSo here are some piece of code I gathered but can't find a way to make them work together
For now my best method is the following based on : http://sgillies.net/blog/870/a-more-perfect-union-continued/. It works fine and I get a list of the 52 states as Shapely geometry. Please feel free to comment if there is a more straight forward way to do this part.
from osgeo import ogrfrom shapely.wkb import loadsfrom numpy import asarrayfrom shapely.ops import cascaded_unionds = ogr.Open('counties.shp')layer = ds.GetLayer(0)#create a list of unique states identifier to be able#to loop through them laterSTATEFP_list = []for i in range(0 , layer.GetFeatureCount()) : feature = layer.GetFeature(i) statefp = feature.GetField('STATEFP') STATEFP_list.append(statefp)STATEFP_list = set(STATEFP_list)#Create a list of merged polygons = states #to be written to filepolygons = []#do the actual dissolving based on STATEFP#and append polygonsfor i in STATEFP_list : county_to_merge = [] layer.SetAttributeFilter("STATEFP = '%s'" %i ) #I am not too sure why "while 1" but it works while 1: f = layer.GetNextFeature() if f is None: break g = f.geometry() county_to_merge.append(loads(g.ExportToWkb())) u = cascaded_union(county_to_merge) polygons.append(u)#And now I am totally stuck, I have no idea how to write #this list of shapely geometry into a shapefile using the#same properties that my source.So the writing is really not straight forward from what I have seen, I really just want the same shapefile only with the country dissolve into states, I don't even need much of the attribute table but I am curious to see how you can pass it on from the source to the new created shapefile.
I found many pieces of code for writing with fiona but I am never able to make it work with my data.Example from How to write Shapely geometries to shapefiles? :
from shapely.geometry import mapping, Polygonimport fiona# Here's an example Shapely geometrypoly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])# Define a polygon feature geometry with one attributeschema = { 'geometry': 'Polygon', 'properties': {'id': 'int'},}# Write a new Shapefilewith fiona.open('my_shp2.shp', 'w', 'ESRI Shapefile', schema) as c: ## If there are multiple geometries, put the "for" loop here c.write({ 'geometry': mapping(poly), 'properties': {'id': 123}, })The problem here is how to do the same with a list of geometry and how to recreate the same properties than the source.Thank you!
أكثر...