Skip to content

Instantly share code, notes, and snippets.

@arthur-e
Created August 7, 2013 17:03
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save arthur-e/6176054 to your computer and use it in GitHub Desktop.
Save arthur-e/6176054 to your computer and use it in GitHub Desktop.

Apparently even a difference of the least precision of the upper left-hand corner causes an alignment error in PostGIS. Consider this example, where I am comparing a land cover raster subset by an MTBS raster and the MTBS raster itself (with the goal being to add the latter to the former, creating the "burned" land cover).

```
SELECT ST_MetaData(r1.rast)
  FROM
(SELECT rast
   FROM geowepp_burnedarea
  WHERE geowepp_burnedarea.rid = 2) r1
UNION
SELECT ST_MetaData(r2.rast)
  FROM
(SELECT ST_Union(ST_Clip(geowepp_landcover.rast, ST_Envelope(geowepp_burnedarea.rast))) AS rast 
   FROM geowepp_landcover, geowepp_burnedarea 
  WHERE ST_Intersects(geowepp_landcover.rast, geowepp_burnedarea.rast)
    AND geowepp_burnedarea.rid = 2) r2
```

And the output of the ST_MetaData()...

```
rid |       upperleftx |       upperlefty | width | height | scalex | scaley | skewx | skewy |  srid | numbands
----+------------------+------------------+-------+--------+--------+--------+-------+-------+-------+----------
1   |    233973.043581 |   4374647.949473 |   482 |    391 |     30 |    -30 |     0 |     0 | 32613 |        1
2   | 233973.043581451 | 4374647.94947295 |   482 |    391 |     30 |    -30 |     0 |     0 | 32613 |        1
```

These two rasters were resampled on the same grid using ArcGIS tools. For practical purposes, they do line up but, as we can see from the output above, they do not line up exactly. At least we don't have to deal with scaling or skew issues! The solution, which is probably best referred to as "fudging it" is to just set the alignment of one of the rasters to that of the other. Here is an example of the solution in action; previously, the ST_MapAlgebraExpr() would have thrown an error because the input rasters were not aligned.

```
SELECT ST_AsTiff(ST_MapAlgebraExpr(ST_SetUpperLeft(r1.rast,
    ST_UpperLeftX(r2.rast), ST_UpperLeftY(r2.rast)), r2.rast, '[rast1.val] + [rast2.val]', NULL, 'UNION')) AS rast
  FROM
(SELECT rast
   FROM geowepp_burnedarea
  WHERE geowepp_burnedarea.rid = 2) r1,
(SELECT ST_Union(ST_Clip(geowepp_landcover.rast, ST_Envelope(geowepp_burnedarea.rast))) AS rast 
   FROM geowepp_landcover, geowepp_burnedarea 
  WHERE ST_Intersects(geowepp_landcover.rast, geowepp_burnedarea.rast)
    AND geowepp_burnedarea.rid = 2) r2
```

Misalignment also causes errors in ST_Union().

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment