Как из точек в CSV получить полигоны в SHP

[Atom74]

Задача в следующем:
Есть трек работы сельхозтехники на поле(точки в WGS-84). У техники есть понятие ширины захвата.
Необходимо получить прямоугольные полигоны, шириной равной ширине захвата, соединяющие эти точки в SHP формате.
Как преобразовать в линии через ogs2ogs и VRT-драйвер разобрался.
Разобрался как рисовать полигоны через WKT.

Как прибавить к точке ширину захвата?
Может есть более "красивый" способ решения данной задачи, например полилиния с указанием её ширины?

[Ariki]

В QGIS Вектор -> Геообработка -> Буферные зоны

[Atom74]

Да, то что надо, только бы в пакетном режиме, без установки и "тыкании" в QGIS.

Ставлю буферную зону равную 0.00005 получаю заливку в 11.1 метров, пополам 5.55 метров.
В чем указывается буферная зона?

Еще проблемка с линиями, необходимо получить не одну линию, а много линий(проходы техники). Признак новой линии есть у точки. Для формирования линии использовал модуль Points2One, но в нем нет такой возможности.
Может есть другой способ формировать линии или другой модуль?

Спасибо.

[trir]

SQL

[Atom74]

Например есть точки:

lat,lon,is_new
10,10,1
20,20,0
30,30,1
40,40,0

Как на SQL получить две линии?

Спасибо.

[trir]

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

Set CrLn = geometryfromtext(concat('LineString(',x(p1),' ',y(p1),',',x(CrP),' ',y(CrP),')'));

Пример создания базы данных
Через OGR загрузить данные в БД (PostGIS, MySQL MS SQL Server), а потом уже работать с ними на SQL

[Atom74]

После преобразования данных к такому виду:
lat,lon,is_new,id_line
10,10,1,1
20,20,0,1
30,30,1,2
40,40,0,2
Модуль Points2One с Group features by id_line нормально строит две линии.
Затем используя совет от Ariki "Вектор -> Геообработка -> Буферные зоны", получаем требуемый результат.

Теперь вопрос возможно ли все это проделать в пакетном режиме без запуска QGIS, например через ogr2ogr?

[trir]

SQL

[Ariki]

Можно с помощью скрипта на Python:
https://pcjericks.github.io/py-gdalogr- ... ate-buffer

[gamm]

Ставлю буферную зону равную 0.00005 получаю заливку в 11.1 метров, пополам 5.55 метров.
В чем указывается буферная зона?

в координатах, в вашем случае - в градусах (0.00005*111000=5.55, 111км о длина градуса на экваторе). А площади у вас будут в квадратных градусах

[trir]

ST_Buffer

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

--A lighter but lamer circle
-- only 2 points per quarter circle is an octagon
--Below is a 100 meter octagon
-- Note coordinates are in NAD 83 long lat which we transform
to Mass state plane meter and then buffer to get measurements in meters;
SELECT ST_AsText(ST_Buffer(
ST_Transform(
ST_SetSRID(ST_MakePoint(-71.063526, 42.35785),4269), 26986),100,2)) As octagon;

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

SELECT ST_Buffer(
 ST_GeomFromText(
  'LINESTRING(50 50,150 150,150 50)'
 ), 10, 'endcap=flat join=round');

[Atom74]

В каком пакете идет библиотека "ogr"?
Использую установщик от OsGeo4W64.
Хотел проверить скрипт https://pcjericks.github.io/py-gdalogr- ... ate-buffer

[Ariki]

Сам не пользуюсь OSGeo4W, но по идее вот это:
pkg-gdal-python

[Atom74]

Сейчас имеется shp файл с точками в атрибутах у которых есть номер принадлежности к линии, необходимо получить shp файл с линиями сгруппированными по номеру точек.

Возможно ли в пакетном режиме использовать модуль Points2One?

Как на Python использовать модуль Points2One?

Спасибо.

[Atom74]

Таким кодом получил из точек линии:

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

import ogr, os

def createLine(inputPointfn, outputLinefn):

# Open shape file with points
    inputds = ogr.Open(inputPointfn)
    inputlyr = inputds.GetLayer()

    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputLinefn):
        shpdriver.DeleteDataSource(outputLinefn)
    outputLineds = shpdriver.CreateDataSource(outputLinefn)

# Create Layer Multiline
    Linelyr = outputLineds.CreateLayer(outputLinefn, geom_type=ogr.wkbMultiLineString)
    featureDefn = Linelyr.GetLayerDefn()
    outFeature = ogr.Feature(featureDefn)

    oMultiLine = ogr.Geometry(ogr.wkbMultiLineString)
    oLine = ogr.Geometry(ogr.wkbLineString)

    for feature in inputlyr:
        oPoint = feature.GetGeometryRef()
	if feature.GetField("field_4") == "1" :
	    oMultiLine.AddGeometry(oLine)
	    oLine = ogr.Geometry(ogr.wkbLineString)
	    oLine.AddPoint(oPoint.GetX(),oPoint.GetY());
	else:
	    oLine.AddPoint(oPoint.GetX(),oPoint.GetY());


    oMultiLine.AddGeometry(oLine)
    outFeature.SetGeometry(oMultiLine)
    Linelyr.CreateFeature(outFeature)


def main(inputPointfn, outputLinefn):
    createLine(inputPointfn, outputLinefn)


if __name__ == "__main__":
    inputPointfn = 'corgon_point.shp'
    outputLinefn = 'corgon_line.shp'

    main(inputPointfn, outputLinefn)

Затем формируем буферные зоны:

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

import ogr, os

def createBuffer(inputfn, outputBufferfn, bufferDist):
    inputds = ogr.Open(inputfn)
    inputlyr = inputds.GetLayer()

    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputBufferfn):
        shpdriver.DeleteDataSource(outputBufferfn)
    outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
    bufferlyr = outputBufferds.CreateLayer(outputBufferfn, geom_type=ogr.wkbPolygon)
    featureDefn = bufferlyr.GetLayerDefn()

    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferDist)
	outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        bufferlyr.CreateFeature(outFeature)


def main(inputfn, outputBufferfn, bufferDist):
    createBuffer(inputfn, outputBufferfn, bufferDist)


if __name__ == "__main__":
    inputfn = 'corgon_line.shp'
    outputBufferfn = 'corgon_poligon.shp'
    bufferDist = 0.000135

    main(inputfn, outputBufferfn, bufferDist)

Как сделать чтобы буферные зоны не сливались в один объект при наложении буферных зон?