Did the earth move for you too?

January 10, 2009 12:24 PM

GISfilter -- I have a PostGIS table containing spatial data in the form of a few thousand shapefiles. Anybody know a quick way I can move *all* of them by about 40 metres?

Due to a bizarre set of circumstances I need to work with a data set of things contructed from GPS in Ireland. Because of a systemic error, every multipolygon I have is out of whack slightly. Boo!

Yet, because the country I live in is small enough for this to work, I can have a bunch of stuff accurate to within 2m if I globally move every multipolygon 50m to the North and 23.4m to the West. Yay!

But I'm too unfamiliar with both Postgres and Python (I'm working with Geodjango for the first time) to figure out which 10-line loop script will do that lateral translation for me. Boo!

So... um... can anybody lend me a hand here? I don't really care which language the one-use-ever script I need is written in. And my multipolygons are strings of co-ords using WGS84, where each X&Y value is separated by a space, and each pair separated by a comma.

The application of your collective genius is warmly appreciated.

posted by genghis to Computers & Internet (11 answers total)

are the numbers being added to X and Y the same? i.e. are you adding 5 to X and 7 to Y, repeated 10,000 times?
posted by mecran01 at 12:33 PM on January 10, 2009

Also, can you post one sample record here?
posted by mecran01 at 12:33 PM on January 10, 2009

Response by poster: a) Yes, that's exactly it, and
b) Not usefully - asking Postgres for one row of the DB shows the multipolygon as one gigantic string of hex values
posted by genghis at 12:40 PM on January 10, 2009

Response by poster: ...and yes, that definitely contradicts my memory of how those values were originally fed to the database. My apologies.
posted by genghis at 12:41 PM on January 10, 2009

In SQL it works to have a query like "UPDATE sometable SET blah = blah + 72" — it will compute the new values for each record based on the previous values. It's been a while since I used the Postgres GIS extensions, but I bet there's a function which takes a polygon and translates it, in which case the answer is a SQL oneliner.

A quick google turns up things like ST_Affine and ST_Translate, but there might be easier ways.
posted by hattifattener at 12:51 PM on January 10, 2009

Best answer: Assuming the geometry is stored as meters north and east of some reference point (probably the Equator and Greenwich), try this:

BEGIN;
UPDATE table SET location = ST_Translate(location, 50, -23.4, 0);
SELECT ST_AsText(location) FROM table;
-- if the locations look OK:
COMMIT;
-- otherwise
ROLLBACK;
-- and try something else

posted by nicwolff at 2:46 PM on January 10, 2009

Response by poster: The values are stored in decimal degrees. Unfortunately, when I tried:

UPDATE table set mpoly = ST_Translate(mpoly, 0.00049, -0.000234, 0);

...the values, while pleasingly altered uniformly, were shifted by a half/a quarter of a degree rather than the much smaller shift we're looking for.
posted by genghis at 4:31 PM on January 10, 2009

Best answer: That's why we did it in a transaction :) Your numbers are a little off though: A degree of longitude is roughly 111km at the equator, so in Ireland at 53ºN it's cos(53/180*π)*111km = 66.8km, and you want to move 23.4m west so that's -.00075º. A degree of latitude is still 111km so 50m/111km = .00045º.

And I got longitude and latitude backward in the query - I think it's x, y, z so let's assume that it's longitude first, so it'd be ST_Translate(location, -.00075, .00045, 0)

No idea why your points are moving so far though - just to test the units, what happens if you ST_Translate(location, 1, 1, 0)?
posted by nicwolff at 9:12 PM on January 10, 2009

Best answer: And what SRID are the polygons defined in? The units of the translate command depend on the SRID, afaict.
posted by hattifattener at 11:33 PM on January 10, 2009

Sorry again! 23.4m / 66.8km = .00035º not .00075º

And yeah, SRID is probably the reason for the scale error.
posted by nicwolff at 12:00 AM on January 11, 2009

Response by poster: SRID=4326, I should have mentioned before.

That more or less worked (it appears the leading zero before the decimal, which I used and you didn't, makes a difference because of a bug).

Y'all rule. AskMefi trinkets now being distributed accordingly. Thanks!
posted by genghis at 1:05 PM on January 11, 2009