Rozpuszczanie wielokątów opartych na atrybutach Pythona (foremna, fiona)?

Próbowałem stworzyć funkcję, która robi w zasadzie to samo, co funkcja „rozpuszczenia” QGIS. Myślałem, że to będzie super łatwe, ale najwyraźniej nie. Więc z tego, co zebrałem, użycie fiona z kształtnym powinno być tutaj najlepszą opcją. Właśnie zacząłem bawić się plikami wektorowymi, więc ten świat jest dla mnie całkiem nowy i python.

Dla tych przykładów pracuję z utworzonym przez hrabstwo plikiem hrabstwa utworzonym tutaj http://tinyurl.com/odfbanu Więc oto fragment kodu, który zebrałem, ale nie mogę znaleźć sposobu, aby mogli współpracować

Na razie moja najlepsza metoda jest następująca: https://sgillies.net/2009/01/27/a-more-perfect-union-continued.html . Działa dobrze i otrzymuję listę 52 stanów jako kształtną geometrię. Prosimy o komentarz, jeśli istnieje prostszy sposób wykonania tej części.

from osgeo import ogr
from shapely.wkb import loads
from numpy import asarray
from shapely.ops import cascaded_union

ds = ogr.Open('counties.shp')
layer = ds.GetLayer(0)

#create a list of unique states identifier to be able
#to loop through them later
STATEFP_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 file
polygons = []

#do the actual dissolving based on STATEFP
#and append polygons
for 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.

Więc pisanie nie jest tak naprawdę proste z tego, co widziałem, naprawdę chcę tylko ten sam plik kształtu tylko po tym, jak kraj rozpłynie się w państwa, nie potrzebuję nawet dużej części tabeli atrybutów, ale jestem ciekawy, jak możesz przejść od źródła do nowo utworzonego pliku kształtu.

Znalazłem wiele fragmentów kodu do pisania za pomocą fiona, ale nigdy nie jestem w stanie sprawić, by działał z moimi danymi. Przykład z Jak pisać geometrie kształtów w plikach kształtów? :

from shapely.geometry import mapping, Polygon
import fiona

# Here's an example Shapely geometry
poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])

# Define a polygon feature geometry with one attribute
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

# Write a new Shapefile
with 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},
    })

Problem polega na tym, jak zrobić to samo z listą geometrii i jak odtworzyć te same właściwości, co źródło.

Pytanie dotyczy Fiony i Shapely, a inna odpowiedź korzystająca z GeoPandów wymaga również znajomości Pand . Ponadto GeoPandas używa Fiony do odczytu / zapisu plików kształtów.

Nie kwestionuję tutaj użyteczności GeoPandas, ale możesz to zrobić bezpośrednio za pomocą Fiony, używając standardowych modułów itertools , szczególnie za pomocą polecenia groupby („W skrócie, groupby bierze iterator i dzieli go na pod-iteratory na podstawie zmian w „kluczu” głównego iteratora. Oczywiście odbywa się to bez wczytywania całego iteratora źródłowego do pamięci, itertools.groupby ).

Oryginalny plik kształtów pokolorowany przez pole STATEFP

from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools
with fiona.open('cb_2013_us_county_20m.shp') as input:
    # preserve the schema of the original shapefile, including the crs
    meta = input.meta
    with fiona.open('dissolve.shp', 'w', **meta) as output:
        # groupby clusters consecutive elements of an iterable which have the same key so you must first sort the features by the 'STATEFP' field
        e = sorted(input, key=lambda k: k['properties']['STATEFP'])
        # group by the 'STATEFP' field 
        for key, group in itertools.groupby(e, key=lambda x:x['properties']['STATEFP']):
            properties, geom = zip(*[(feature['properties'],shape(feature['geometry'])) for feature in group])
            # write the feature, computing the unary_union of the elements in the group with the properties of the first element in the group
            output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})

Wynik

Bardzo polecam GeoPandas do radzenia sobie z dużym asortymentem funkcji i wykonywania operacji masowych.

Rozszerza ramki danych Pandas i używa zgrabnie pod maską.

from geopandas import GeoSeries, GeoDataFrame

# define your directories and file names
dir_input = '/path/to/your/file/'
name_in = 'cb_2013_us_county_20m.shp'
dir_output = '/path/to/your/file/'
name_out = 'states.shp'

# create a dictionary
states = {}
# open your file with geopandas
counties = GeoDataFrame.from_file(dir_input + name_in)

for i in range(len(counties)):
    state_id = counties.at[i, 'STATEFP']
    county_geometry = counties.at[i, 'geometry']
    # if the feature's state doesn't yet exist, create it and assign a list
    if state_id not in states:
        states[state_id] = []
    # append the feature to the list of features
    states[state_id].append(county_geometry)

# create a geopandas geodataframe, with columns for state and geometry
states_dissolved = GeoDataFrame(columns=['state', 'geometry'], crs=counties.crs)

# iterate your dictionary
for state, county_list in states.items():
    # create a geoseries from the list of features
    geometry = GeoSeries(county_list)
    # use unary_union to join them, thus returning polygon or multi-polygon
    geometry = geometry.unary_union
    # set your state and geometry values
    states_dissolved.set_value(state, 'state', state)
    states_dissolved.set_value(state, 'geometry', geometry)

# save to file
states_dissolved.to_file(dir_output + name_out, driver="ESRI Shapefile")

Jako dodatek do @ genu odpowiedzi , musiałem rozpuścić przez więcej niż jednego pola tak zmodyfikował swój kod do obsługi wielu pól. Poniższy kod wykorzystuje operator.itemgetterdo grupowania według wielu pól:

# Modified from /gis//a/150001/2856
from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools
from operator import itemgetter


def dissolve(input, output, fields):
    with fiona.open(input) as input:
        with fiona.open(output, 'w', **input.meta) as output:
            grouper = itemgetter(*fields)
            key = lambda k: grouper(k['properties'])
            for k, group in itertools.groupby(sorted(input, key=key), key):
                properties, geom = zip(*[(feature['properties'], shape(feature['geometry'])) for feature in group])
                output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})


if __name__ == '__main__':
    dissolve('input.shp', 'input_dissolved.shp', ["FIELD1", "FIELD2", "FIELDN"))