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

33

Digamos que armé un Shapefile y que todas las características tienen sus vértices cambiados 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 cambio 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 la función agradable ShiftCoords . Sin embargo, el hilo fue muy informativo! ¡Gracias a todos!

    
pregunta Jose 05.03.2012 - 17:42

10 respuestas

21

Utilizando 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";
    
respondido por el David Bitner 05.03.2012 - 20:40
28

He diseñado Fiona (un contenedor de OGR) para hacer este tipo de procesamiento simple.

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

    
respondido por el sgillies 05.03.2012 - 19:13
13

Y aunque la publicación está etiquetada con python, dado 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);
});
    
respondido por el Tim Schaub 05.03.2012 - 21:01
11

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.

    
respondido por el Antonio Falciano 02.10.2013 - 18:05
8

El GRASS GIS v.edit module :

Se supone una ubicación y un conjunto de mapas existentes 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)
    
respondido por el webrian 08.03.2012 - 08:09
8

Otra opción sería utilizar las opciones de reproyección simplemente en ogr2ogr, sin duda un enfoque más intrépido que los enfoques de JEQL, Fiona o GeoScript, pero de todos modos eficaz. Tenga en cuenta que las proyecciones desde y hasta no tienen que ser realmente la proyección real del shapefile original, siempre que lo único que esté cambiando entre las proyecciones utilizadas en s_srs y t_srs sean el falso 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 eliminó 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 ejecuta:

ogr2ogr --optfile projcs.txt shifted.shp input.shp
    
respondido por el David Bitner 06.03.2012 - 19:16
7

Una opción R que usa el paquete maptools y su función de 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 = ""))
}
    
respondido por el mdsumner 07.03.2012 - 01:19
4

Usando el analizador de shapefile en geofunciones, podría usar XSLT para realizar el proceso. Por supuesto, necesitarías volver a convertir el archivo de forma 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>
    
respondido por el Peter Rushforth 07.03.2012 - 21:56
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
    })
}  
    
respondido por el jericks 10.03.2012 - 23:06
0

Aquí está la versión de OGR

driver = ogr.GetDriverByName ("ESRI Shapefile")

def mover (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   
    
respondido por el Moshe Yaniv 27.09.2017 - 07:55

Lea otras preguntas en las etiquetas