Bagaimana cara mudah menggeser semua fitur dalam dataset vektor?

33

Katakanlah saya mengumpulkan Shapefile dan semua fitur memiliki simpulnya yang digeser dengan jumlah yang konstan. Apa cara termudah untuk menggeser semua fitur (karenanya posisi (x, y) dari simpulnya) oleh pergeseran sewenang-wenang? Saya memiliki banyak file yang akan saya gunakan koreksi ini, jadi jawaban Bash / OGR akan lebih disukai :)

Akhirnya, saya akhirnya menggunakan Spatialite untuk ini, karena memiliki fungsi yang bagus ShiftCoords. Namun, utas itu sangat informatif! Terima kasih semuanya!

Jose
sumber
Saya suka entri ini. Seluruh halaman ini adalah contoh luar biasa dari T&J yang berjalan dengan benar. Sebuah pertanyaan yang dijelaskan dengan jelas, dan setiap jawaban memberikan solusi yang unik, valid, dan lengkap. Cantiknya. Saya mengubah setiap suara, masing-masing berdasarkan kemampuan mereka sendiri.
matt wilkie
@Jose Posting ini perlu sedikit pembaruan karena peningkatan yang relatif baru dari perpustakaan GDAL. Ada solusi satu liner sekarang (lihat jawaban di bawah)! Dimungkinkan untuk menggunakan fungsi SpatiaLite ShiftCoords secara langsung dengan utilitas ogr2ogr.
Antonio Falciano

Jawaban:

21

Menggunakan JEQL Ini dapat dilakukan dengan tiga baris:

ShapefileReader t file: "shapefile.shp";
out = select * except (GEOMETRY), Geom.translate(GEOMETRY,100,100) from t;
ShapefileWriter out file: "ahapefile_shift.shp";
David Bitner
sumber
canggih! bagus!
WolfOdrade
Bagus, David. Itu ketat.
sgillies
1
Hanya harus menunjukkan ... "ahapefile?"
WolfOdrade
Saya akhirnya menggunakan fungsi terjemahan Spatialite, yang mirip dengan apa yang Anda sarankan di sini.
Jose
30

Saya telah merancang Fiona (pembungkus OGR) untuk membuat pemrosesan seperti ini sederhana.

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)

Pembaruan : Saya telah meletakkan versi yang berbeda, lebih kencang dari skrip ini di http://sgillies.net/blog/1128/geoprocessing-for-hipsters-translating-features .

sgillies
sumber
2
"Geoprosesing untuk para hipsters" Saya berharap saya bisa menaikkan peringkat 10x itu untuk gelar yang luar biasa itu
Ragi Yaser Burhum
13

Dan meskipun postingannya ditandai dengan python, karena JEQL telah disebutkan, inilah contoh dengan JavaScript (menggunakan 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
sumber
13

Menggunakan GDAL> = 1.10.0 dikompilasi dengan SQLite dan SpatiaLite:

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

di mana shiftX = 1 dan shiftY = 10.

Antonio Falciano
sumber
1
Brilliant - solusi CLI satu baris sederhana.
Dave X
pendek dan mudah!
Kurt
8

Modul GRASS GIS v.edit :

Lokasi dan mapset yang ada dalam proyeksi pencocokan diasumsikan.

Dalam skrip shell:

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

atau dalam skrip 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
sumber
8

Opsi lain adalah menggunakan opsi proyeksi ulang hanya di ogr2ogr, tentu saja pendekatan peretasan daripada pendekatan JEQL, Fiona, atau GeoScript tetapi efektif. Perhatikan bahwa proyeksi dari dan ke tidak benar-benar harus menjadi proyeksi aktual dari bentuk asli selama satu-satunya hal yang berubah antara proyeksi yang digunakan dalam s_srs dan t_srs adalah arah timur dan utara palsu. Dalam contoh ini saya hanya menggunakan Google Mercator. Saya yakin ada sistem koordinat yang jauh lebih sederhana untuk digunakan sebagai basis, tapi yang ini tepat di depan saya untuk menyalin / menempel.

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

Atau untuk menyimpan pengetikan / tempel, simpan yang berikut ke projcs.txt(sama seperti di atas, tetapi dihapus dengan menyertakan tanda kutip tunggal):

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

lalu jalankan:

ogr2ogr --optfile projcs.txt shifted.shp input.shp
David Bitner
sumber
2
Ini berubah menjadi golf geo-scripting! Langkah selanjutnya adalah meretas tabel EPSG Anda untuk menghilangkan PROJCS literal yang panjang;)
sgillies
@sgillies, tidak perlu meretas epsg, cukup simpan projcs ke file dan gunakan --optfile, mis ogr2ogr --optfile projcs.txt shifted.shp input.shp. Saya akan melipatnya menjadi jawaban.
matt wilkie
7

Opsi R menggunakan paket maptools dan fungsi elide-nya:

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

Menggunakan parser shapefile di fungsi-fungsi geografis, Anda bisa menggunakan XSLT untuk melakukan proses. Tentu saja Anda perlu mengonversi kembali ke shapefile sesudahnya :-).

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

Ini adalah versi GeoScript Groovy:

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
    })
}  
tersentak
sumber
0

Ini adalah versi OGR

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
sumber