These are my notes on understanding/evaluating raster support in Postgis 2.0 .
Platform
- Windows 7
- PostgreSQL9.0
- Postgis2.0
Key Resources
- Installation on Windows:
Jorge Arevalo’s notes are key. (But expect ‘minor’ variations based on your own system’s configurations, esp. with gdal) - Loading rasters into postgis tables and getting some understanding of it
- Pierre Racine’s WKT Raster tutorial (Replace gdal2wktraster.py with raster2pgsql.py)
- Jorge Arevalo’s notes
- This example from Fuzzy Tolerance
- Working with postgis raster: intersections with vector datasets
1. Example data i used
- SRTM for southern India downloaded from here, clipped in QGIS using gdaltools plugin, called srtm_clip.tif
- A vector layer i had of administrative boundaries, called wards
2. Loading an elevation raster into Postgis2.0
a. Using raster2pgsql.py, in a command line, like so:
raster2pgsql.py -r srtm_clip.tif -s 4326 -I -t srtm_table -k 700×700 -o srtm_clip.sql
where
- srtm_clip.tif is the input elevation raster
- 4326 is the EPSG code for lat-long, wgs 84 coordinate system
- -I tells it to make a spatial index
- -k 700×700 is the size of each tile to break the input raster into
- srtm_sql is the output file that will contain the sql commands that creates the raster table in the database that has previously been set up
b.Load the sql file
Use pgAdminIII’s sql editor to open the ‘srtm_clip.sql’ the above step created, and run it.
FYI this is what the generated sql file looked like. (and i could not get the psql command line option for loading the sql file to work). It ran in less than a second.
BEGIN;
CREATE TABLE “public”.”srtm_table” (rid serial PRIMARY KEY);
SELECT AddRasterColumn(‘public’,’srtm_table’,’rast’,4326, ARRAY[’16BSI’], false, true, null, 0.000833333333333, -0.000833333333333, 700, 700, ST_Envelope(ST_SetSRID(‘POLYGON((77.256195534800000 13.681944121300001,77.256195534800000 12.561944121300002,78.131195534800000 13.681944121300001,78.131195534800000 12.561944121300002,77.256195534800000 13.681944121300001))’::geometry, 4326)));
INSERT INTO “public”.”srtm_table” ( rast ) VALUES ( (‘01000001004F1BE8….’)::raster );
INSERT INTO “public”.”srtm_table” ( rast ) VALUES ( (‘01000001004F1BE8….’)::raster );
INSERT INTO “public”.”srtm_table” ( rast ) VALUES ( (‘01000001004F1BE8….’)::raster );
INSERT INTO “public”.”srtm_table” ( rast ) VALUES ( (‘01000001004F1BE8….’)::raster );
CREATE INDEX “srtm_table_rast_gist_idx” ON “public”.”srtm_table” USING GIST (st_convexhull(rast));
END;
c. What do you get?
- Refresh the database and you’ll see a new table called srtm_table.
- It has two columns rid and rast, and 4 rows. each row is for a rid, which i guess is a tile id.
- Each rid references one tile – choosing k=700×700 caused the input raster to be split into 4 tiles
- Smaller tiles are generally better – but i did run into one issue (you’ll see in just a bit)
Viewing the tiles (in OpenJump, because Qgis 1.6 does not work with Postgis2.0, even for non-rasters).
4 srtm tiles (overlaid over administrative boundaries)
-
Useful commands at this point
— postgis versions
select postgis_full_version(), postgis_raster_lib_build_date(), postgis_raster_lib_version();— raster formats
SELECT short_name, long_name
FROM st_gdaldrivers();–metadata
SELECT (md).*, (bmd).*
FROM (SELECT ST_Metadata(rast) AS md,
ST_BandMetadata(rast) AS bmd
FROM srtm_table LIMIT 1
) foo;
— no data value was not set even though i did put it into the QGIS-gdal plugin application–size
SELECT pg_size_pretty(pg_total_relation_size(‘srtm_table’));
–gives 1.8mb
—
–no-data value setting
UPDATE srtm_table SET rast = ST_SetBandNoDataValue(rast,-32768);
vacuum analyze srtm_table;–summary stats (mean) of the raster
SELECT ST_Quantile(rast, 0.5) As value
FROM srtm_table;
— this seems to give me the median of each rid (but i want the median of the entire raster)
1;886
2;339
3;857
4;0
2. Intersecting raster with vector layers
a. performing the intersection
This below worked, creating 70,503 rows – but it took 80 minutes!!!
create table srtm_wards_inter as
select
ward_name,ward_id, (gv).geom As the_geom,(gv).val
FROM (Select ward_name, ward_id,ST_Intersection(rast,geom) As gv
From srtm_table, wardscleanattributes
WHERE geom && rast
and ST_Intersects(rast, geom)) foo;
b. extracting area-weighted mean elevations per ward polygon (sorted by mean elevation)
SELECT ward_id,ward_name,
sum(ST_Area(ST_Transform(the_geom, 32643)) * val) /
sum(ST_Area(ST_Transform(the_geom, 32643))) AS meanelev
FROM srtm_wards_inter
GROUP BY ward_id, ward_name
ORDER BY meanelev, ward_id;
–gives me the ward id’s, names and mean elevations for all 198 wards
159;”Kengeri”;792.638912308088
131;”Nayandahalli”;821.362355078139
130;”Ullalu”;828.12323070928
198;”Hemmigepura”;831.628739711782
…
67;”Nagapura”;921.610375368031
46;”Jayachamarajendra Nagar”;922.392379372865
179;”Shakambarinagar”;922.91608088907
62;”Ramaswamy Palya”;923.071218071726
9;”Vidyaranyapura”;923.304108544661
16;”Jalahalli”;924.921270389213
63;”Jayamahal”;926.957157345059
93;”Vasanthnagar”;928.315490103474
194;”Gottigere”;931.443435808278
34;”Gangenahalli”;932.719869810253
45;”Malleshwaram”;932.814742613641
35;”Aramanenagar”;935.117666378835
Major Questions
- Why did the intersection take so long (see 2a above)
- How can i get raster statistics irrespective of rid?
- example: quantile stats for the raster are split by rid (i.e tile)
- example: if i get median elevation by ward like so:
- SELECT ward_id,ward_name,ST_Quantile(rast,0.5) As median
- FROM srtm_table,wardscleanattributes
- WHERE ST_Intersects(rast,geom)
- ORDER BY ward_id, median;
I STILL get multiple entries per ward, if (i) the ward crosses two or more rid (tiles), or (ii) the ward is split into 2 polygons
If anyone can help with the above, I can update my post soon.
Thanks!
Vishal Mehta's weblog 