Центроиды линий

[Игорь Лебедь]

Как получить центры линий (точки, лежащие на линии и равноудалённые от её краёв)? Центроиды полигонов работают только для полигональных слоёв. Не понятно, работает ли функция centroid для линий? И работает ли он в шейп-файлах или же только в базе?

[Александр Мурый]

Можно попробовать через питоновский модуль shapely, как написано здесь. Или через GRASS (там есть пару дополн. модулей для получения средних точек).

[Игорь Лебедь]

Можно попробовать через питоновский модуль shapely, как написано здесь. Или через GRASS (там есть пару дополн. модулей для получения средних точек).

1. А как их установить, и нет ли GUI-аналога?
2. Какие конкретно модули GRASS, не подскажете?
Ещё питоновский модуль считает только для одной линии, как изменить его для всех линий слоя? И как вывести результат не в принт, а в колонку? Пробовал так:

Код: Выделить всё

"xypoint"=(midPoint.to_wkt())

не получилось

[Александр Мурый]

Рекомендуемый модуль для GRASS 7 - v.centerpoint, ставится в GUI через установщик дополнений (см. в меню Settings - Addons extensions).

[Игорь Лебедь]

Спасибо. Но что-то я его не нашёл в списке, есть v.line.center. А в QGIS в панель инструментов GRASS нельзя его установить?

Спойлер

Безымянный.png (46.04 КБ) 11509 просмотров

[Александр Мурый]

Спасибо. Но что-то я его не нашёл в списке, есть v.line.center. А в QGIS в панель инструментов GRASS нельзя его установить?

У вас точно GRASS 7 стоит? Потому как <v.line.center> работает только под GRASS 6.4.* (и только под Linux, т.к. это шелл-скрипт), а <v.centerpoint> - только под GRASS 7.

[Александр Мурый]

Вот простой Питон-скрипт с использованием библиотек <OGR> и <shapely>, на входе - линейный шейп-файл, на выходе - точечный с "серединными" точками.

Код: Выделить всё

center_points.py lines.shp cent_pts.shp

Код: Выделить всё

#!/usr/bin/env python

# -*- coding: utf-8 -*-



import sys, os

from osgeo import ogr

from shapely import wkt

from shapely.geometry import LineString



in_lns = sys.argv[1]

out_pts = sys.argv[2]



# open shapefile with OGR

input = ogr.Open(in_lns)

in_lyr = input.GetLayer(0)



## create output point shapefile

ogrDriver = ogr.GetDriverByName("ESRI Shapefile")

outDataSource = ogrDriver.CreateDataSource(out_pts)



# copy SRS from the input shapefile

srs = in_lyr.GetSpatialRef()



# create layer

out_lyr = outDataSource.CreateLayer(out_pts, srs, ogr.wkbPoint)



# export to WKT, compute center points and write them to shp

for f in in_lyr:

    g = f.geometry().ExportToWkt()

    shapelyLine = LineString(wkt.loads(g))

    midPoint = shapelyLine.interpolate(shapelyLine.length/2).to_wkt()

    # write point features to layer

    feature = ogr.Feature(out_lyr.GetLayerDefn())

    point = ogr.CreateGeometryFromWkt(midPoint)

    feature.SetGeometry(point)

    out_lyr.CreateFeature(feature)

    feature.Destroy()



# Destroy the data source to free resources

outDataSource.Destroy() 

center_points.py

(1.05 КБ) 541 скачивание

[Игорь Лебедь]

Спасибо. Но что-то я его не нашёл в списке, есть v.line.center. А в QGIS в панель инструментов GRASS нельзя его установить?

У вас точно GRASS 7 стоит? Потому как <v.line.center> работает только под GRASS 6.4.* (и только под Linux, т.к. это шелл-скрипт), а <v.centerpoint> - только под GRASS 7.

Тут вот 7, а десктопная 6,4:

Спойлер

Безымянный2.png (21.76 КБ) 11472 просмотра

меня интересует-то главным образом эта, панелька. Она зависит от десктопной или нет?

[Игорь Лебедь]

Большое спасибо за скрипт, но не могу понять, как запустить его? Можно пример с конкретным синтаксисом путей в Windows? И в какой консоли его запускать? Консоль питон в кугис/коммандная строка в осгеов/cmd?

Вот простой Питон-скрипт с использованием библиотек <OGR> и <shapely>, на входе - линейный шейп-файл, на выходе - точечный с "серединными" точками.

Код: Выделить всё

center_points.py lines.shp cent_pts.shp

Результат:
C:\OSGeo4W>center_points.py StreetsKK.shp points.shp
Fatal Python error: Py_Initialize: unable to load the file system codec
File "C:\OSGeo4W\\apps\Python27\lib\encodings\__init__.py", line 123
raise CodecRegistryError,\
^
SyntaxError: invalid syntax

[Александр Мурый]

меня интересует-то главным образом эта, панелька. Она зависит от десктопной или нет?

С взаимоотношениями QGIS и GRASS 7 без пол-литра не разберешься. Как-то добавлять модули в Processing можно, но как именно, быстро не подскажу. Сам такого пока не делал.

Можно пример с конкретным синтаксисом путей в Windows? И в какой консоли его запускать? Консоль питон в кугис/коммандная строка в осгеов/cmd?

Скрипт запускается в командной строке (в Windows через cmd.exe). Для его работы должны быть установлены: GDAL и питоновские привязки к нему (см. здесь для Win), а также библиотека Shapely (см., например, здесь). Все это для Питона 2.7, не для 3.*

Переменная PATH для пользователя должна включать GDAL+Python+Python scripts
Например, для WinXP у меня так

Код: Выделить всё

C:\Program Files\GDAL;C:\Python27;C:\Python27\Scripts 

P.S. для Windows желательно убрать из скрипта две первые строки (начинающиеся с символа #).

[Игорь Лебедь]

Установил всё как указали, запускаю из cmd, ругается:

Код: Выделить всё

No module named osgeo

Как понимаю, не правильно задана переменная PATH, но в ней прописан путь к установленному в OSGEO Python (на Windows 7 64):

Код: Выделить всё

C:\Program Files\GDAL;C:\OSGeo4W64\apps\Python27\;C:\OSGeo4W64\apps\Python27\Scripts\; ***

Как определить папку, которую нужно прописать в переменную?
При попытке запустить из панели инструментов вообще ругается:

Спойлер

qwsd.jpg (94.33 КБ) 11333 просмотра

[Александр Мурый]

Насчёт OSGeo4W и панели инструментов не подскажу. Замечу только, что у вас в PATH намешаны отдельно установленный GDAL с OSGeo4W - вряд ли это хорошо.

[Александр Мурый]

Вот более простой и логичный вариант без использования Shapely. Открываем файл в редакторе Питон-консоли или делаем пользовательский скрипт на панели инструментов; ну или вставляем нижеприведённый код куда надо; выделяем в списке слоёв слой с линиями; запускаем скрипт. На выходе - точечный memory-файл.

Код: Выделить всё

from qgis.core import *

import qgis.utils



layer = qgis.utils.iface.mapCanvas().currentLayer()

iter = layer.getFeatures()



crs = layer.crs().toWkt()



out_name = layer.name() + '_centerpoints'

vl = QgsVectorLayer("Point?crs=" + crs, out_name, "memory")

pr = vl.dataProvider()



fet = QgsFeature()



for feature in iter:

    geom = feature.geometry()

    length = geom.length()

    mid = geom.interpolate(length/2)

    fet.setGeometry(mid)

    pr.addFeatures([fet])



QgsMapLayerRegistry.instance().addMapLayer(vl)

center_points_pyqgis.py

(518 байт) 479 скачиваний

[Игорь Лебедь]

Спасибо! Заработало! Однако, в этом случае в слое не сохраняются атрибуты исходного линейного, что желательно получить, и связи с исходным слоем тоже не работают. Можно ли как-то таким методом унаследовать родительские атрибуты?
При попытке запуска первого скрипта в редакторе консоли питон выходит ошибка:

Код: Выделить всё

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "G:/center_points.py", line 6, in <module>
    in_lns = sys.argv[1]
IndexError: list index out of range

[Александр Мурый]

Однако, в этом случае в слое не сохраняются атрибуты исходного линейного, что желательно получить, и связи с исходным слоем тоже не работают. Можно ли как-то таким методом унаследовать родительские атрибуты?

Как-то можно, естественно. Надо разбираться. Просто это не входило в изначальную задачу.

При попытке запуска первого скрипта в редакторе консоли питон выходит ошибка:

Это нормально, т.к. первый скрипт не предназначен для запуска в Питон-консоли, только в системном терминале.