Область пересечения двух шейп-файлов - Python

У меня есть два шейп-файла на Python, и я хотел бы найти площадь всех пространств, которые они перекрывают.

Я могу использовать sjoin из геопанд, чтобы получить области, в которых они соединяются, но для мест, где есть несколько перекрытий, я хотел бы оставить только одно с наибольшей площадью.

municipality = gpd.read_file(muni_file)
soil_type = gpp.read_file(soil)
combined = gpd.sjoin(municipality,soil_type,how="left",op="intersects")

С помощью OGR я могу получить площадь многоугольника, как показано ниже.

from osgeo import ogr

wkt = "POLYGON ((1162440.5712740074 672081.4332727483, 1162440.5712740074 647105.5431482664, 1195279.2416228633 647105.5431482664, 1195279.2416228633 672081.4332727483, 1162440.5712740074 672081.4332727483))"
poly = ogr.CreateGeometryFromWkt(wkt)

Поэтому мне интересно, есть ли способ взять мой объединенный шейп-файл и получить область, где они пересекаются, чтобы я оставил только максимальное значение для каждого муниципалитета.


python geopandas
person DarknessFalls    schedule 01.06.2017    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Да, я считаю, что вы можете объединить цикл через применить и получить размер каждого из пересечений.

начните с переиндексации комбинированного, поскольку я предполагаю, что они дублируются из sjoin ()

combined = combined.reset_index()

Затем определите вспомогательную функцию (get_size_of_intersection), затем мы объединим цикл, применим get_size_of_intersection () и создадим новую серию с именем crossction_size.

некоторые примечания:

-комбинированный будет иметь геометрию муниципалитета

-combined будет иметь столбец / ряд с именем index_right, который будет индексом почвы_типа.

-поскольку мы имеем дело с красивыми объектами, мы можем использовать перекресток () и свойство area.

def get_size_of_intersection(row, soil_type):
    return row['geometry'].intersection(soil_type['geometry'].iloc[int(row['index_right'])]).area

combined['intersection_size'] = combined.apply(lambda row : 
                                       get_size_of_intersection(row, soil_type), axis=1)

Мы создадим еще одну серию под названием max_intersection_size. Здесь я предполагаю, что у муниципалитета есть какая-то серия 'name', по которой мы можем сгруппировать и применить max ()

combined['max_intersection_size'] = combined.groupby('name')['intersection_size'].transform(max)

Затем, используя логическую индексацию, мы получаем желаемые данные

(т.е. где размер_пересечения равен max_intersection_size)

filter = combined['intersection_size'] == combined['max_intersection_size']
combined[filter]
person Bob Haffner    schedule 10.06.2017