¿Cómo cambiar fácilmente todas las funciones en un conjunto de datos vectoriales?

33

Digamos que armé un Shapefile y todas las entidades tienen sus vértices desplazados en una cantidad constante. ¿Cuál es la forma más fácil de cambiar todas las características (de ahí la posición (x, y) de sus vértices) mediante un desplazamiento arbitrario? Tengo muchos archivos a los que aplicaría esta corrección, por lo que sería preferible una respuesta Bash / OGR :)

Finalmente, terminé usando Spatialite para esto, ya que tiene una buena función ShiftCoords. Sin embargo, el hilo fue muy informativo! ¡Gracias a todos!

Jose
fuente
Me encanta esta entrada. Toda esta página es un excelente ejemplo de preguntas y respuestas correctas. Una pregunta clara y clara, y cada respuesta ofrece una solución única, válida y completa. Es bonito. Voté a cada uno de ellos, cada uno por sus propios méritos.
Matt Wilkie
@Jose Esta publicación necesita una pequeña actualización debido a las mejoras relativamente recientes de la biblioteca GDAL. ¡Ahora hay una solución de una línea (vea la respuesta a continuación)! Es posible utilizar la función SpatiaLite ShiftCoords directamente con la utilidad ogr2ogr.
Antonio Falciano

Respuestas:

21

Usando JEQL Esto se puede hacer con tres líneas:

ShapefileReader t file: "shapefile.shp";
out = select * except (GEOMETRY), Geom.translate(GEOMETRY,100,100) from t;
ShapefileWriter out file: "ahapefile_shift.shp";
David Bitner
fuente
¡innovador! ¡bonito!
WolfOdrade
Bien David. Está apretado.
sgillies
1
Solo tengo que señalar ... "ahapefile?"
WolfOdrade
Terminé usando la función de traducción de Spatialite, que es similar a lo que sugieres aquí.
Jose
30

Diseñé Fiona (un contenedor OGR) para simplificar este tipo de procesamiento.

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)

Actualización : he puesto una versión diferente y más ajustada de este script en http://sgillies.net/blog/1128/geoprocessing-for-hipsters-translating-features .

sgillies
fuente
2
"Geoprocesamiento para los hipsters" Ojalá pudiera votar 10 veces por ese título increíble
Ragi Yaser Burhum,
13

Y aunque la publicación está etiquetada con python, ya que JEQL ya se ha mencionado, aquí hay un ejemplo con JavaScript (usando 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
fuente
13

Usando GDAL> = 1.10.0 compilado con SQLite y SpatiaLite:

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

donde shiftX = 1 y shiftY = 10.

Antonio Falciano
fuente
1
Brillante: una solución simple de CLI de una línea.
Dave X
corto y fácil!
Kurt
8

El módulo GRASS GIS v.edit :

Se supone una ubicación y un conjunto de mapas existentes en la proyección coincidente.

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)
Webrian
fuente
8

Otra opción sería usar las opciones de reproyección simplemente en ogr2ogr, ciertamente un enfoque más hackear que los enfoques JEQL, Fiona o GeoScript, pero no obstante efectivos. Tenga en cuenta que las proyecciones desde y hacia realmente no necesitan ser la proyección real del archivo de forma original siempre que lo único que esté cambiando entre las proyecciones utilizadas en s_srs y t_srs sean la falsa este y el norte. En este ejemplo, solo estoy usando Google Mercator. Estoy seguro de que hay un sistema de coordenadas mucho más simple para usar como base, pero este estaba justo frente a mí para copiar / pegar.

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

O para guardar la escritura / pegado, guarde lo siguiente en projcs.txt(igual que arriba, pero quitado entre comillas simples):

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

y luego ejecuta:

ogr2ogr --optfile projcs.txt shifted.shp input.shp
David Bitner
fuente
2
¡Se está convirtiendo en golf de geo-script! El siguiente paso sería hackear su tabla EPSG para eliminar los PROJCS literales largos;)
sgillies
@sgillies, no es necesario hackear epsg, solo guarde los proyectos en un archivo y úselos --optfile, por ejemplo ogr2ogr --optfile projcs.txt shifted.shp input.shp. Lo doblaré en la respuesta.
Matt Wilkie
7

Una opción R usando 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 = ""))
}
mdsumner
fuente
4

Usando el analizador shapefile en geofunciones, puede usar XSLT para realizar el proceso. Por supuesto, necesitaría volver a convertirlo en 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>
Peter Rushforth
fuente
4

Aquí hay una versión 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
    })
}  
imbéciles
fuente
0

Aquí está la versión de OGR

driver = ogr.GetDriverByName ("ESRI Shapefile")

movimiento def (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
fuente