Version 6 (modified by 12 years ago) ( diff ) | ,
---|
In morphology, the dilation operation followed by erosion is an extremely useful technique that allows for the simplification of complex shapes. For example, the technique can be used to produce a simplified coastline, where bays and fjords are ignored, and that can later be used as a mask for further operations. Below is an illustration of the process of this technique for simplifying the convoluted coastline of Newfoundland, Canada, using GSHHS.
At first glance, it would seem that simply using one of the lower resolution coastline datasets from GSHHS would achieve the same result, but they were generated using a different procedure, which is not equivalent.
Import GSHHS polygons into PostGIS
GSHHS shapefiles provided by the maintainers can be imported into a PostGIS database (named gshhs here) using:
shp2pgsql -s 4326 -I -g geom GSHHS_f_L1.shp gshhs_f_l1 | psql -d gshhs
Subset coastlines
Before performing any dilation/erosion, it is helpful to subset the coastline polygons for efficiency purposes. The gid
field in the table holding the GSHHS data identifies a single polygon. The island of Newfoundland has gid=23
:
CREATE OR REPLACE VIEW newfoundland_subset_f SELECT gshhs_f_l1.gid, gshhs_f_l1.geom FROM gshhs_f_l1 WHERE gshhs_f_l1.gid = 23;
Note that more subtle subsetting can be done by interesecting with user-defined polygons and other geometries to obtain the exact desired polygon.
Dilation and erosion
PostGIS can perform these two operations via ST_Buffer
. For the purposes outlined here, the process involves buffering positively (dilating) the full resolution coastline by a given amount k
, and then buffering the result negatively (eroding) by the same amount. The following shell script applies the technique using different values for k
. It uses the utility ogr2ogr
to export the PostGIS objects to an ASCII table in GMT (Generic Mapping Tools) format:
COAST="newfoundland_subset_f.gmt" BUFOUT="newfoundland_subset_f_bufout" BUFIN="newfoundland_subset_f_bufin" BUFFERS=( 0.08 0.10 0.12 0.14 0.16 ) . gmt_shell_functions.sh # from GMT distribution TMPDIR=/var/tmp gmt_init_tmpdir ogr2ogr -f GMT ${COAST} \ PG:"host=localhost user=username password=password dbname=gshhs" \ ${COAST%%.gmt} for buf in ${BUFFERS[@]}; do cat <<EOF > ${GMT_TMPDIR}/buffers.sql CREATE OR REPLACE VIEW ${COAST%%.gmt}_bufout AS SELECT gid, ST_Buffer(geom, ${buf}) FROM ${COAST%%.gmt}; CREATE OR REPLACE VIEW ${COAST%%.gmt}_bufin AS SELECT gid, ST_Buffer(st_buffer, -${buf}) FROM ${COAST%%.gmt}_bufout; EOF psql -f ${GMT_TMPDIR}/buffers.sql gshhs ogr2ogr -f GMT ${COAST%%.gmt}_bufout${buf}.gmt \ PG:"host=localhost user=username password=password dbname=gshhs" \ ${COAST%%.gmt}_bufout ogr2ogr -f GMT ${COAST%%.gmt}_bufin${buf}.gmt \ PG:"host=localhost user=username password=password dbname=gshhs" \ ${COAST%%.gmt}_bufin done # Clean up gmt_remove_tmpdir
The animated GIF
image below illustrates the results:
Attachments (2)
- newfoundland.2.gif (116.1 KB ) - added by 12 years ago.
- newfoundland.gif (180.6 KB ) - added by 12 years ago.
Download all attachments as: .zip