Como mudar facilmente todos os recursos em um conjunto de dados vetoriais?

33

Digamos que eu montei um Shapefile e todos os recursos tenham seus vértices alterados em uma quantidade constante. Qual é a maneira mais fácil de mudar todas as características (daí a posição (x, y) de seus vértices) por uma mudança arbitrária? Eu tenho muitos arquivos aos quais eu aplicaria essa correção, portanto, uma resposta Bash / OGR seria preferida :)

Por fim, acabei usando o Spatialite para isso, pois tem a boa função ShiftCoords. No entanto, o tópico foi muito informativo! Obrigado a todos!

Jose
fonte
Eu simplesmente amo esta entrada. Esta página inteira é um excelente exemplo de perguntas e respostas corretas. Uma pergunta direta e claramente descrita, e cada resposta fornece uma solução única, válida e completa. É lindo. Votei cada um deles, cada um por seus próprios méritos.
Matt Wilkie
@ Jose Este post precisa de uma pequena atualização devido às melhorias relativamente recentes da biblioteca GDAL. Agora existe uma solução única (veja a resposta abaixo)! É possível usar a função SpatiaLite ShiftCoords diretamente com o utilitário ogr2ogr.
Antonio Falciano

Respostas:

21

Usando JEQL Isso pode ser feito com três linhas:

ShapefileReader t file: "shapefile.shp";
out = select * except (GEOMETRY), Geom.translate(GEOMETRY,100,100) from t;
ShapefileWriter out file: "ahapefile_shift.shp";
David Bitner
fonte
ponta! legais!
21412 WolfOdrade
Legal David. Isso é apertado.
sgillies
1
Só tenho que apontar ... "ahapefile?"
7892 WolfOdrade
Acabei usando a função de tradução do Spatialite, que é semelhante ao que você sugere aqui.
318 Jose Jose
30

Eu projetei o Fiona (um wrapper OGR) para simplificar esse tipo de processamento.

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)

Atualização : coloquei uma versão diferente e mais rígida desse script em http://sgillies.net/blog/1128/geoprocessing-for-hipsters-translating-features .

sgillies
fonte
2
"Geoprocessamento para descolados" Eu gostaria de poder apenas aprovar esse 10x por esse título incrível
Ragi Yaser Burhum 08/03/12
13

E apesar de a postagem ter sido marcada com python, já que o JEQL já foi mencionado, aqui está um exemplo com 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
fonte
13

Usando GDAL> = 1.10.0 compilado com SQLite e SpatiaLite:

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

onde shiftX = 1 e shiftY = 10.

Antonio Falciano
fonte
1
Brilhante - uma solução CLI de uma linha simples.
Dave X
curto e fácil!
Kurt
8

O módulo GRASS GIS v.edit :

Um local e um conjunto de mapas existentes na projeção correspondente são assumidos.

Em um 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

ou em um script 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
fonte
8

Outra opção seria usar as opções de reprojeção simplesmente no ogr2ogr, certamente uma abordagem mais obscena do que as abordagens JEQL, Fiona ou GeoScript, mas eficazes, no entanto. Observe que as projeções from e to realmente não precisam ser a projeção real do shapefile original, desde que a única coisa que esteja mudando entre as projeções usadas nos s_srs e t_srs sejam o leste e o norte falsos. Neste exemplo, estou apenas usando o Google Mercator. Estou certo de que existe um sistema de coordenadas muito mais simples para usar como base, mas este estava bem na minha frente para copiar / colar.

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

Ou, para salvar a digitação / colagem, salve o seguinte em projcs.txt(o mesmo que acima, mas removido aspas 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"]]

e então execute:

ogr2ogr --optfile projcs.txt shifted.shp input.shp
David Bitner
fonte
2
Está se transformando em golfe de geo-script! O passo seguinte seria a de cortar sua mesa EPSG para eliminar os PROJCS longo literais;)
sgillies
@sgillies, não há necessidade de hackear o epsg, basta salvar os projetos em um arquivo e usar --optfile, por exemplo ogr2ogr --optfile projcs.txt shifted.shp input.shp. Vou dobrá-lo na resposta.
Matt Wilkie
7

Uma opção R usando pacotes maptools e sua função 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
fonte
4

Usando o analisador shapefile em geofunções, você pode usar o XSLT para executar o processo. Claro que você precisaria converter novamente em shapefile depois :-).

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

Aqui está uma versão do 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
    })
}  
idiotas
fonte
0

Aqui está a versão OGR

driver = ogr.GetDriverByName ("ESRI Shapefile")

def movimento (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
fonte