Desplazamiento de todas las características en el conjunto de datos vectoriales mediante Bash/OGR

Question

Supongamos que he creado un archivo Shapefile y todas las características tienen sus vértices desplazados por una cantidad constante.

¿Cuál es la forma más fácil de desplazar todas las características (de ahí el (x,y) posición de sus vértices) por un desplazamiento arbitrario?

Tengo muchos archivos a los que aplicaría esta corrección, así que sería preferible una respuesta de Bash/OGR.

---

Terminé usando Spatialite para esto, ya que tiene la agradable función ShiftCoords .

Answer 1

El SIG de GRASS v.editar módulo :

Se supone que existe una ubicación y un conjunto de mapas en la proyección correspondiente.

En un script de 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

o en un script de 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)

Answer 2

Una opción de R que utiliza el paquete maptools y su función 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 = ""))
}

Answer 3

Utilizando el analizador de shapefile en geofunciones, podría utilizar XSLT para realizar el proceso. Por supuesto, tendría que volver a convertir a shapefile después :-).

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

Answer 4

Aquí hay una versión de 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
    })
}

Answer 5

Aquí está la versión 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