Clipping Satellite Images at the Coastline

Raster data from satellite born sensors will, in most cases, cover the land areas and oceans as well. If the satellite data is intended only for background imagery, then having the sea areas appear is fine. However, when the data will be put thru some GIS analysis, the sea areas can actually interfere.  Here’s a simple method, using GRASS GIS to clip out only the land surface.

Two sources of global elevation data are now available online. The well known SRTM project offers a DEM covering much of the globe, at a resolution of 90 meters, and includes bathymetric depths as well. More recently, the ASTER sensor aboard NASA’s  Earth Observing Satellite has produced ASTER GDEM. This dataset covers even more of the earth’s surface than SRTM and at a resolution of 30 meters. The ASTER data also includes values beyond the shorelines. However since this sensor worked differently than SRTM, the sea floor depth is not necessarily correct. If the DEM is used to delineate streams and watersheds, then anomalies can appear in coastal regions.   The best way to use this DEM is to change all sea surface areas to a NULL value.  We can use NOAA’s Global shoreline (GSHHS) data to clip out the sea surface areas with the following steps.

For this exercise, I’ll get an ASTER tile along the coastline of Croatia, with plenty of islands. For those who have never used ASTER data, the procedure for getting data is somewhat involved. You must be registered on their WIST site, and logged in. Then a multi-step search process guides you to the data tiles you want. Eventually, after committing your request, you’ll get an email with an ftp link for downloading the zipped files. So, get registered on the WIST site and search for and download the ASTER DEM tile for your area.

Now download the “Global Self-consistent, Hierarchical, High-resolution Shoreline” shapefile from NOAA. Here’s what the original ASTER DEM looks like in Quantum GIS, overlayed with the GSHHS layer:

Croatia coastline

The GSHHS data are split into separate shapefiles at different resolutions. The highest resolution is designated ‘f’ (full). The next highest appears as ‘h’ (high). Once you unzip the downloaded file, the shapefile components appear in subdirectories named with those same single letters. So the highest resolution shoreline is in the ‘f’ directory.

Now we import the data into GRASS. Start GRASS in a LOCATION defined as EPSG:4326.  First import the ASTER DEM (and bne sure to set the current region to match that raster):
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia > \ in=~/geodata/data/ASTGTM_N43E016_dem.tif out=aster_dem
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia> g.region -p rast=aster_dem
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia> r.colors aster_dem color=srtm
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia> r.shaded.relief aster_dem shade=aster_shade

Next we want to bring in the GSHHS shapefile. This is a pretty large dataset, and we’ll only need a very small portion of the coastline so we utilize an option to the GRASS module “spatial=” to select right from the start only those parts of the GSHHS data that intersect our region of interest.
(Note that there is a GRASS Addon module specifically for importing the GSHHS data in their original netCDF binary format, and projecting to the current LOCATION on the way. We’re demonstrating here by simply importing the GHSSH shapefile.)

You should be in a WGS84, EPSG:4326 LOCATION. I choose the spatial parameters to just cover the extent of the ASTER tile. We get the values from the output of the above g.region -p line. The commands look like this:
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia > \ dsn=~/geodata/data/GSHHS_shp/f/GSHHS_f_L1.shp out=gshhs spatial="15.9,42.9,17.1,44.1"
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia> gshhs
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia> v.centroids gshhs out=gshhs_areas
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia> gshhs_areas out=MASK use=val value=1
GRASS 6.4.0 (WGS84):~/geodata/grass/WGS84/Croatia> r.mapcalc aster_dem_land=aster_dem

What we’ve done is convert the GSHHS areas (land areas) to a raster named MASK, then using r.mapcalc a copy of the ASTER DEM is made which, due to the mask, will have the original DEM values on land surfaces, and NULL over sea surfaces.  And here’s the result:

aster clipped
Resulting DEM with shaded relief

Leave a Reply

Your email address will not be published. Required fields are marked *