Wie können alle Features in einem Vektordatensatz einfach verschoben werden?

33

Angenommen, ich habe ein Shapefile zusammengestellt und bei allen Features sind die Scheitelpunkte um einen konstanten Betrag verschoben. Was ist der einfachste Weg, um alle Merkmale (daher die (x, y) -Position ihrer Scheitelpunkte) durch eine willkürliche Verschiebung zu verschieben? Ich habe viele Dateien, auf die ich diese Korrektur anwenden würde, daher wäre eine Bash / OGR-Antwort vorzuziehen :)

Schließlich habe ich Spatialite dafür verwendet, da es die nette Funktion hat ShiftCoords. Der Thread war jedoch sehr informativ! Vielen Dank an alle!

Jose
quelle
Ich liebe diesen Eintrag. Diese ganze Seite ist ein hervorragendes Beispiel für eine fehlerfreie F & A. Eine einfach zu beschreibende Frage und jede Antwort ergibt eine eindeutige, gültige und vollständige Lösung. Es ist wunderschön. Ich habe jeden einzelnen für sich gestimmt.
Matt Wilkie
@Jose Dieser Beitrag benötigt ein kleines Update, da die GDAL-Bibliothek in relativ kurzer Zeit erweitert wurde. Es gibt jetzt eine einzeilige Lösung (siehe Antwort unten)! Es ist möglich, die SpatiaLite ShiftCoords-Funktion direkt mit dem Dienstprogramm ogr2ogr zu verwenden.
Antonio Falciano

Antworten:

21

Verwenden von JEQL Dies kann in drei Zeilen erfolgen:

ShapefileReader t file: "shapefile.shp";
out = select * except (GEOMETRY), Geom.translate(GEOMETRY,100,100) from t;
ShapefileWriter out file: "ahapefile_shift.shp";
David Bitner
quelle
innovativ, auf dem neuesten Stand! nett!
WolfOdrade
Schön, David. Das ist eng.
sgillies
1
Ich muss nur darauf hinweisen ... "ahapefile?"
WolfOdrade
Am Ende habe ich die Übersetzungsfunktion von Spatialite verwendet, die Ihren Vorschlägen hier ähnelt.
Jose
30

Ich habe Fiona (einen OGR-Wrapper) entwickelt, um diese Art der Verarbeitung zu vereinfachen.

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)

Update : Ich habe eine andere, engere Version dieses Skripts unter http://sgillies.net/blog/1128/geoprocessing-for-hipsters-translating-features veröffentlicht .

sgillies
quelle
2
"Geoverarbeitung für Hipster" Ich wünschte, ich könnte das 10-fache für diesen tollen Titel
verbessern
13

Und obwohl der Beitrag mit Python getaggt ist, ist hier ein Beispiel mit JavaScript (unter Verwendung von GeoScript ) , da JEQL bereits erwähnt wurde .

/**
 * 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
quelle
13

Mit GDAL> = 1.10.0 kompiliert mit SQLite und SpatiaLite:

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

Dabei ist shiftX = 1 und shiftY = 10.

Antonio Falciano
quelle
1
Brilliant - eine einfache CLI-Lösung mit einer Leitung.
Dave X
kurz und einfach!
Kurt
8

Das Modul GRASS GIS v.edit :

Ein vorhandener Ort und eine vorhandene Karte in der übereinstimmenden Projektion werden angenommen.

In einem Shell-Skript:

#!/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

oder in einem Python-Skript:

#!/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
quelle
8

Eine andere Möglichkeit wäre, die Reprojektionsoptionen einfach in ogr2ogr zu verwenden, was sicherlich hackiger ist als die JEQL-, Fiona- oder GeoScript-Ansätze, aber trotzdem effektiv. Beachten Sie, dass die Projektionen von und bis nicht unbedingt die eigentliche Projektion des ursprünglichen Shapefiles sein müssen, solange sich zwischen den in s_srs und t_srs verwendeten Projektionen nur die falsche Ost- und Nordrichtung ändert. In diesem Beispiel verwende ich nur Google Mercator. Ich bin mir sicher, dass es ein viel einfacheres Koordinatensystem gibt, das ich als Basis verwenden kann, aber dieses war genau vor mir, um es zu kopieren / einzufügen.

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

Oder speichern Sie Folgendes, um das Eingeben / Einfügen zu speichern projcs.txt(wie oben, jedoch ohne einfache Anführungszeichen):

-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"]]

und dann laufen:

ogr2ogr --optfile projcs.txt shifted.shp input.shp
David Bitner
quelle
2
Es verwandelt sich in Geo-Scripting-Golf! Der nächste Schritt wäre, Ihre EPSG-Tabelle zu hacken, um die lange wörtliche PROJCS zu beseitigen;)
sgillies
@sgillies, keine Notwendigkeit, EPSG zu hacken, speichern Sie einfach die Projekte in eine Datei und verwenden Sie --optfile, z ogr2ogr --optfile projcs.txt shifted.shp input.shp. Ich werde es in die Antwort falten.
Matt Wilkie
7

Eine R-Option unter Verwendung des Pakets maptools und seiner elide-Funktion:

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
quelle
4

Mit dem Shapefile-Parser in Geofunktionen können Sie den Prozess mithilfe von XSLT ausführen. Natürlich müsstest du danach wieder in Shapefile konvertieren :-).

<?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
quelle
4

Hier ist eine Groovy GeoScript-Version:

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
quelle
0

Hier ist die OGR-Version

driver = ogr.GetDriverByName ("ESRI-Shapefile")

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   
Moshe Yaniv
quelle