wiki:UsersWikiBufferCoast

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"
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}.gmt \
    PG:"host=localhost user=sluque password=otariidae dbname=gshhs" \
    ${COAST}
for buf in ${BUFFERS[@]}; do
    cat <<EOF > ${GMT_TMPDIR}/buffers.sql
    CREATE OR REPLACE VIEW ${COAST}_bufout AS
    SELECT gid, ST_Buffer(geom, ${buf}) FROM ${COAST};
    CREATE OR REPLACE VIEW ${COAST}_bufin AS
    SELECT gid, ST_Buffer(st_buffer, -${buf}) FROM ${COAST}_bufout;
EOF
    psql -f ${GMT_TMPDIR}/buffers.sql gshhs
    ogr2ogr -f GMT ${COAST}_bufout${buf}.gmt \
	PG:"host=localhost user=sluque password=otariidae dbname=gshhs" \
	${COAST}_bufout
    ogr2ogr -f GMT ${COAST}_bufin${buf}.gmt \
	PG:"host=localhost user=sluque password=otariidae dbname=gshhs" \
	${COAST}_bufin
done

# Clean up
gmt_remove_tmpdir

The animated GIF image below illustrates the results:

Last modified 5 years ago Last modified on Nov 27, 2012 7:15:01 PM

Attachments (2)

Download all attachments as: .zip