Jak łatwo przesunąć wszystkie funkcje w zbiorze danych wektorowych?

33

Powiedzmy, że utworzyłem plik kształtu i wszystkie elementy mają przesunięte wierzchołki o stałą wartość. Jaki jest najprostszy sposób na przesunięcie wszystkich cech (stąd pozycja (x, y) ich wierzchołków) o dowolne przesunięcie? Mam wiele plików, do których chciałbym zastosować tę poprawkę, więc preferowana byłaby odpowiedź Bash / OGR :)

W końcu użyłem do tego Spatialite, ponieważ ma on fajną funkcję ShiftCoords. Wątek był jednak bardzo pouczający! Dziękuje wszystkim!

Jose
źródło
Po prostu uwielbiam ten wpis. Ta cała strona jest doskonałym przykładem odpowiedzi na pytania i odpowiedzi. Proste, jasno opisane pytanie, a każda odpowiedź daje unikalne, ważne i kompletne rozwiązanie. To jest piękne. Głosowałem za każdym z nich osobno.
matt wilkie
@Jose Ten post wymaga niewielkiej aktualizacji ze względu na stosunkowo nowe ulepszenia biblioteki GDAL. Jest teraz jedno rozwiązanie liniowe (patrz odpowiedź poniżej)! Możliwe jest użycie funkcji SpatiaLite ShiftCoords bezpośrednio za pomocą narzędzia ogr2ogr.
Antonio Falciano

Odpowiedzi:

21

Korzystanie z JEQL Można to zrobić za pomocą trzech linii:

ShapefileReader t file: "shapefile.shp";
out = select * except (GEOMETRY), Geom.translate(GEOMETRY,100,100) from t;
ShapefileWriter out file: "ahapefile_shift.shp";
David Bitner
źródło
ostrze! miły!
WolfOdrade
Fajnie, David. To jest ciasne.
sgillies,
1
Muszę tylko wskazać ... „ahapefile?”
WolfOdrade,
Skończyło się na użyciu funkcji translacji Spatialite, która jest podobna do tego, co tutaj sugerujesz.
Jose
30

Zaprojektowałem Fionę (owijkę OGR), aby ten rodzaj przetwarzania był prosty.

from fiona import collection
import logging

log = logging.getLogger()

# A few functions to shift coords. They call eachother semi-recursively.
def shiftCoords_Point(coords, delta):
    # delta is a (delta_x, delta_y [, delta_y]) tuple
    return tuple(c + d for c, d in zip(coords, delta))

def shiftCoords_LineString(coords, delta):
    return list(shiftCoords_Point(pt_coords, delta) for pt_coords in coords)

def shiftCoords_Polygon(coords, delta):
    return list(
        shiftCoords_LineString(ring_coords, delta) for ring_coords in coords)

# We'll use a map of these functions in the processing code below.
shifters = {
    'Point': shiftCoords_Point,
    'LineString': shiftCoords_LineString,
    'Polygon': shiftCoords_Polygon }

# Example 2D shift, 1 unit eastward and northward
delta = (1.0, 1.0)

with collection("original.shp", "r") as source:

    # Create a sink for processed features with the same format and 
    # coordinate reference system as the source.
    with collection(
            "shifted.shp", 
            "w",
            driver=source.driver,
            schema=source.schema,
            crs=source.crs
            ) as sink:

        for rec in source:
            try:
                g = rec['geometry']
                g['coordinates'] = shifters[g['type']](
                    g['coordinates'], delta )
                rec['geometry'] = g
                sink.write(rec)
            except Exception, e:
                log.exception("Error processing record %s:", rec)

Aktualizacja : Umieściłem inną, ściślejszą wersję tego skryptu na stronie http://sgillies.net/blog/1128/geoprocessing-for-hipsters-translating-features .

sgillies
źródło
2
„Geoprocessing dla hipstersów” Chciałbym po prostu głosować 10 razy za ten niesamowity tytuł
Ragi Yaser Burhum
13

I chociaż post jest otagowany pytonem, ponieważ JEQL został już wspomniany, oto przykład z JavaScript (przy użyciu GeoScript ).

/**
 * Shift all coords in all features for all layers in some directory
 */

var Directory = require("geoscript/workspace").Directory;
var Layer = require("geoscript/layer").Layer;

// offset for all geometry coords
var dx = dy = 10;

var dir = Directory("./data");
dir.names.forEach(function(name) {
    var orig = dir.get(name);
    var shifted = Layer({
        schema: orig.schema.clone({name: name + "-shifted"})
    });
    orig.features.forEach(function(feature) {
        var clone = feature.clone();
        clone.geometry = feature.geometry.transform({dx: dx, dy: dy});
        shifted.add(clone);
    });
    dir.add(shifted);
});
Tim Schaub
źródło
13

Używając GDAL> = 1.10.0 skompilowanych z SQLite i SpatiaLite:

ogr2ogr data_shifted.shp data.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,1,10) FROM data"

gdzie shiftX = 1 i shiftY = 10.

Antonio Falciano
źródło
1
Brilliant - proste jedno liniowe rozwiązanie CLI.
Dave X
krótkie i łatwe!
Kurt
8

Moduł GRASS GIS v.edit :

Zakładana jest istniejąca lokalizacja i zestaw map w pasującej projekcji.

W skrypcie powłoki:

#!/bin/bash

for file in `ls | grep \.shp$ | sed 's/\.shp$//g'`
do
    v.in.ogr dsn=./${file}.shp output=$file
    v.edit map=$file tool=move move=1,1 where="1=1"
    v.out.ogr input=$file type=point,line,boundary,area dsn=./${file}_edit.shp
done

lub w skrypcie Python:

#!/usr/bin/env python

import os
from grass.script import core as grass

for file in os.listdir("."):
    if file.endswith(".shp"):
        f = file.replace(".shp","")
        grass.run_command("v.in.ogr", dsn=file, output=f)
        grass.run_command("v.edit", map=f, tool="move", move="1,1", where="1=1")
        grass.run_command("v.out.ogr", input=f, type="point,line,boundary,area", dsn="./%s_moved.shp" % f)
Webrian
źródło
8

Inną opcją byłoby użycie opcji reprojection po prostu w ogrod2ogr, z pewnością bardziej skomplikowane podejście niż w przypadku JEQL, Fiona lub GeoScript, ale skuteczne mimo wszystko. Zauważ, że rzuty do i do tak naprawdę nie muszą być rzeczywistym rzutem oryginalnego pliku kształtu, o ile jedyną rzeczą, która zmienia się między rzutami używanymi w s_srs i t_srs, jest fałszywe wschodzenie i północ. W tym przykładzie używam tylko Google Mercator. Jestem pewien, że istnieje znacznie prostszy układ współrzędnych, który można wykorzystać jako bazę, ale ten był przede mną do skopiowania / wklejenia.

ogr2ogr -s_srs EPSG:900913 -t_srs 'PROJCS["Google Mercator",GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["semi_minor",6378137.0],PARAMETER["latitude_of_origin",0.0],PARAMETER["central_meridian",0.0],PARAMETER["scale_factor",1.0],PARAMETER["false_easting",1000.0],PARAMETER["false_northing",1000.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","900913"]]' -f "ESRI Shapefile" shift.shp original.shp

Lub, aby zapisać pisanie / wklejanie, zapisz następujące elementy projcs.txt(tak samo jak powyżej, ale usunąłem pojedyncze cudzysłowy):

-s_srs EPSG:900913 -t_srs PROJCS["Google Mercator",GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["semi_minor",6378137.0],PARAMETER["latitude_of_origin",0.0],PARAMETER["central_meridian",0.0],PARAMETER["scale_factor",1.0],PARAMETER["false_easting",1000.0],PARAMETER["false_northing",1000.0],UNIT["m",1.0],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","900913"]]

a następnie uruchom:

ogr2ogr --optfile projcs.txt shifted.shp input.shp
David Bitner
źródło
2
Zmienia się w golfa geo-skryptowego! Kolejnym krokiem byłoby zhakowanie stołu EPSG, aby wyeliminować długi dosłowny PROJCS;)
sgillies
@ sgillies, nie trzeba hakować epsg, po prostu zapisz projcs do pliku i użyj --optfilenp ogr2ogr --optfile projcs.txt shifted.shp input.shp. Złożę to do odpowiedzi.
matt wilkie
7

Opcja R wykorzystująca pakiety maptools i jej funkcję elide:

shift.xy <- c(1, 2)
library(maptools)
files <- list.files(pattern = "shp$")
for (fi in files) {
  xx <- readShapeSpatial(fi)
  ## update the geometry with elide arguments
  shifted <- elide(xx, shift = shift.xy)
  ## write out a new shapfile
  writeSpatialShape(shifted, paste("shifted", fi, sep = ""))
}
mdsumner
źródło
4

Używając parsera plików shapefile w geofunkcjach, możesz użyć XSLT do wykonania tego procesu. Oczywiście musisz później przekonwertować plik shapefile :-).

<?xml version="1.0" encoding="UTF-8"?>
<xsl:stylesheet xmlns:xsl="http://www.w3.org/1999/XSL/Transform"
    version="2.0" xmlns:gml="http://www.opengis.net/gml">
    <xsl:param name="x_shift" select="0.0"/>
    <xsl:param name="y_shift" select="0.0"/>

    <!-- first the identity transform makes sure everything gets copied -->
    <xsl:template match="node()|@*">
        <xsl:copy>
            <xsl:apply-templates select="@*|node()"/>
        </xsl:copy>
    </xsl:template>
    <!-- for any element with coordinate strings, apply the translation factors -->
    <!-- note that a schema-aware processor could use the schema type names to simplify -->
    <xsl:template match="gml:pos|gml:posList|gml:lowerCorner|gml:upperCorner">
        <xsl:copy>
            <!-- this xpath parses the ordinates, assuming x y ordering (shapefiles), applies translation factors -->
            <xsl:value-of select="
                for $i in tokenize(.,'\s+') return 
                  if ($i[(position() mod 2) ne 0]) then 
                    number($i)+$x_shift 
                  else 
                    number($i)+$y_shift
             "/>
        </xsl:copy>
    </xsl:template>
</xsl:stylesheet>
Peter Rushforth
źródło
4

Oto wersja Groovy GeoScript:

import geoscript.workspace.Directory
import geoscript.layer.Layer

int dx = 10
int dy = 10

def dir = new Directory("./data")
dir.layers.each{name ->
    def orig = dir.get(name)
    def shifted = dir.create("${name}-shifted", orig.schema.fields)
    shifted.add(orig.cursor.collect{f ->
        f.geom = f.geom.translate(dx, dy)
        f
    })
}  
jericks
źródło
0

Oto wersja OGR

driver = ogr.GetDriverByName („Plik kształtu ESRI”)

def move (dx, dy, dz):

dataSource = driver.Open(path,1)
layer = dataSource.GetLayer(0)
for feature in layer:
    get_poly = feature.GetGeometryRef()
    get_ring = get_poly.GetGeometryRef(0)
    points   = get_ring.GetPointCount()
    set_ring = ogr.Geometry(ogr.wkbLinearRing)
    for p in xrange(points):
        x,y,z = get_ring.GetPoint(p)
        x += dx
        y += dy
        z += dz
        set_ring.AddPoint(x,y)
        print x,y,z
set_poly = ogr.Geometry(ogr.wkbPolygon)
set_poly.AddGeometry(set_ring)

feature.SetGeometry(set_poly)
layer.SetFeature(feature)

del layer
del feature
del dataSource   
Mosze Yaniv
źródło