Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Let's say that I put together a Shapefile and all the features have their vertices shifted by a constant amount. What's the easiest way of shifting all the features (hence the (x,y) position of their vertices) by an arbitrary shift? I have lots of files that I would apply this correction to, so a Bash/OGR answer would be preferred :)

Finally, I ended up using Spatialite for this, as it has the nice function ShiftCoords. However, the thread was very informative! Thanks all!

share|improve this question
    
I just love this entry. This whole page is a superb example of Q&A gone right. A straightforward clearly described question, and every answer gives a unique, valid and complete solution. It's beautiful. I upvoted every single one, each on their own merits. –  matt wilkie Dec 4 '12 at 6:03
    
@Jose This post needs a little update due to the relatively recent enhancements of GDAL library. There's a one liner solution now (see the answer below)! It's possible to use the SpatiaLite ShiftCoords function directly with the ogr2ogr utility. –  afalciano Oct 3 '13 at 10:18

9 Answers 9

up vote 18 down vote accepted

Using JEQL This can be done with three lines:

ShapefileReader t file: "shapefile.shp";
out = select * except (GEOMETRY), Geom.translate(GEOMETRY,100,100) from t;
ShapefileWriter out file: "ahapefile_shift.shp";
share|improve this answer
    
cutting edge! nice! –  WolfOdrade Mar 5 '12 at 19:45
    
Nice, David. That's tight. –  sgillies Mar 6 '12 at 16:30
1  
Just have to point out..."ahapefile?" –  WolfOdrade Mar 7 '12 at 3:45
    
I ended up using Spatialite's translate function, which is similar to what you suggest here. –  Jose Mar 8 '12 at 17:40

I've designed Fiona (an OGR wrapper) to make this kind of processing 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)

Update: I've put a different, tighter version of this script at http://sgillies.net/blog/1128/geoprocessing-for-hipsters-translating-features.

share|improve this answer
1  
"Geoprocessing for hipsters" I wish I could just upvote that 10x for that awesome title –  Ragi Yaser Burhum Mar 8 '12 at 1:33

And though the post is tagged with python, since JEQL has already been mentioned, here's an example with JavaScript (using 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);
});
share|improve this answer

Another option would be to use the reprojection options simply in ogr2ogr, certainly a hackier approach than either the JEQL, Fiona, or GeoScript approaches but effective none-the-less. Note that the from and to projections don't really need to be the actual projection of the original shapefile as long as the only thing that is changing between the projections used in the s_srs and t_srs are the false easting and northing. In this example I am just using the Google Mercator. I am sure there is a much simpler coordinate system to use as the base, but this one was right in front of me to copy/paste.

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

Or to save typing/pasting, save the following to projcs.txt (same as above, but removed enclosing single quotes):

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

and then run:

ogr2ogr --optfile projcs.txt shifted.shp input.shp
share|improve this answer
2  
It's turning into geo-scripting golf! Next step would be to hack your EPSG table to eliminate the long literal PROJCS ;) –  sgillies Mar 6 '12 at 18:52
    
@sgillies, no need to hack epsg, just save the projcs to a file and use --optfile, e.g. ogr2ogr --optfile projcs.txt shifted.shp input.shp. I'll fold it into the answer. –  matt wilkie Dec 4 '12 at 5:48

An R option using package maptools and its elide function:

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 = ""))
}
share|improve this answer

The GRASS GIS v.edit module:

An existing location and mapset in the matching projection is assumed.

In a shell script:

#!/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

or in a Python script:

#!/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)
share|improve this answer

Using the shapefile parser in geofunctions, you could use XSLT to perform the process. Of course you would need to convert back to shapefile afterwards :-).

<?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>
share|improve this answer

Here is a Groovy GeoScript version:

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
    })
}  
share|improve this answer

Using GDAL >= 1.10.0 compiled with SQLite and SpatiaLite:

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

where shiftX = 1 and shiftY = 10.

share|improve this answer
    
Brilliant -- a simple one line CLI solution. –  Dave X yesterday

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.