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

He diseñado Fiona (una envoltura de 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 .

Answer 2

Utilizando JEQL Esto puede hacerse 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";

Answer 3

Utilizando 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.

Answer 4

Y aunque el post está etiquetado con python, ya que se ha mencionado JEQL, 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);
});

Answer 5

Otra opción sería utilizar las opciones de reproyección simplemente en ogr2ogr, sin duda un enfoque más complicado que los enfoques JEQL, Fiona o GeoScript, pero eficaz en cualquier caso. Tenga en cuenta que las proyecciones desde y hasta no necesitan ser la proyección real del shapefile original, siempre y cuando lo único que cambie entre las proyecciones utilizadas en el s_srs y el t_srs sean la falsa easting y northing. En este ejemplo sólo estoy utilizando el Google Mercator. Estoy seguro de que hay un sistema de coordenadas mucho más simple para utilizar como base, pero este estaba justo en frente de 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 ahorrarse el teclear/pegar, guarde lo siguiente en projcs.txt (igual que el anterior, pero eliminando las 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 correr:

ogr2ogr --optfile projcs.txt shifted.shp input.shp