Coastline extraction using Landsat-8 multispectral imagery and band ratio technique

Prerequisites

• This sample demonstrates how the coastlines can be extracted using ArcGIS Image Server. Alternatively, this can be done using ArcGIS Image for ArcGIS Online.

Introduction

Due to anthropogenic activities and natural processes i.e. changes in sea level, sedimentation and wave energy coastlines are changing throughout the world. Coastline is contact line between sea and land, it is an important linear feature on earth's surface with dynamic nature.

Traditionally, the coastlines were manually digitized which was time-consuming and labour intensive. Remote sensing is a good alternative to extract coastlines using satellite imagery, this way both temporal and spatial aspects can be covered.

Satellite imagery of visible range can be used for interpretation and can be easily obtained. But the imageries covering infrared wavelength is best to extract boundary between land and water. So, the satellites which covers both visible and infrared spectrum are widely accepted for coastline extraction and mapping.

Landsat-8 multispectral imagery is used in the current study as it covers a wavelength ranging from 0.43 to 12.51 micrometers, and hence suitable for coastal and aerosol studies.

Neccessary Imports

import os
import glob
from zipfile import *

import arcgis
import arcpy
from arcpy.management import PolygonToLine
from datetime import datetime
import pandas as pd
from arcgis.features import GeoAccessor, GeoSeriesAccessor
from arcgis.raster.analytics import convert_raster_to_feature
from arcgis.features.manage_data import overlay_layers
from arcgis.raster.functions import equal_to, greater_than, greater_than_equal, clip, apply, extract_band, stretch

Connect to your GIS

from arcgis import GIS
gis =  GIS('home')
gis2 = GIS(profile="your_enterprise_profile")

Get the data for analysis

Multispectral Landsat includes Landsat GLS and Landsat 8 imagery for use in visualization and analysis. This layer is time enabled and includes a number band combinations and indices rendered on demand. The Landsat 8 imagery includes eight multispectral bands from the Operational Land Imager (OLI) with 30m spatial resolution and two bands from the Thermal Infrared Sensor (TIRS) of 100m spatial resolution. It is updated daily with new imagery directly sourced from the Landsat on AWS collection.

landsat_item = gis.content.get('d9b466d6a9e647ce8d1dd5fe12eb434b')
landsat = landsat_item.layers[0]
landsat_item

Multispectral Landsat
Landsat multispectral and multitemporal imagery with on-the-fly renderings and indices for visualization and analysis. The Landsat 8 and 9 imagery in this layer is updated daily and is directly sourced from the USGS Landsat collection on AWS.
Imagery Layer by esri
Last Modified: July 01, 2022
3 comments, 1336187 views

The buffer feature of 20 km was created from USA boundaries. This feature layer geometry will be used to get the Landsat-8 tiles of coastal areas.

aoi = gis.content.get('6ee0006608704d0389881e0bcf936778')
aoi

usa_coast_polygon
usa_coast_polygon
Feature Layer Collection by api_data_owner
Last Modified: January 24, 2022
0 comments, 737 views

Prepare data for analysis

A map widget was created, which will define the extent of area of interest for the analysis.

m = gis.map('USA')
m.content.add(aoi)
m

m.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
 'xmin': -10124000.92647267,
 'ymin': 3341135.188124092,
 'xmax': -9823144.78314236,
 'ymax': 3463434.4333803146}

Multispectral imagery layers consist of data for whole world. First step, is to filter out the cloud free data for the study area which will be used in the analysis.

Create geometry of aoi

The aoi has four polygons representing: Islands, East coast, West coast and Alaska.

dfm = aoi.layers[0].query(out_fields="OBJECTID, type, code").sdf
dfm

In the current analysis the area of interest is situated at East coast. The geometry of aoi was created for filtering out the landsat-8 tiles for the study area. objectid=2 was used as it represents the east coast.

aoi_layer = aoi.layers[0]

## use the correct objectid by refering to the above dataframe for your area of interest. 
aoi_feature = aoi_layer.query(where="objectid=2")
aoi_geom = aoi_feature.features[0].geometry
aoi_geom['spatialReference'] = {'wkid':3857}

Filter out the Landsat-8 tiles

The landsat-8 tiles were filtered out on the basis of date of acquisition and cloud cover. The order of tiles is on the basis of cloud cover from low to high.

import pandas as pd
from datetime import datetime
selected = landsat.filter_by(where="(Category = 1) and (CloudCover<0.03)",
                             time=[datetime(2016, 1, 1), datetime(2019, 12, 31)],
                             geometry=arcgis.geometry.filters.intersects(aoi_geom))

df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear", 
                    order_by_fields="CloudCover").sdf
df['AcquisitionDate'] = pd.to_datetime(df['AcquisitionDate'], unit='ms')
df

## Create a list of tiles 
oid = df["OBJECTID"].tolist()

Create mosaic raster from tiles

After creating a list of landsat-8 tiles, next step was to create a mosaic raster for the study area. mosaic_by was used to lock all the tiles of landsat-8 list.

landsat.mosaic_by(method='lock_rasters', lock_rasters=oid)
landsat.extent = m.extent
landsat

Preprocessing of the data

Extract bands from Landsat-8 tiles

From the mosaic raster, 3 single band raster were created for NIR, Red and Blue band using extract_band function.

Landsat-8 band spectrum table

Landsat-8 fifth band covers 0.85 - 0.88 micrometers of EMR spectrum. NIR band is best for delineating land and water interface. Due to high absorption of water in NIR spectrum the reflectance of water is around zero whereas vegetation and land has high reflectance value. In NIR imagery the water appears black and land & vegetation appears bright grey to white which makes it easier to delineate land and water boundary.

nir1 = extract_band(landsat, [5])
nir1

Second band of landsat-8 represents blue spectrum and covers 0.45 - 0.51 micrometer. Water has high reflectance in blue spectrum.

blue1 = extract_band(landsat, [2])
blue1

Fourth band represents red spectrum and covers 0.64 - 0.67 micrometer wavelength of EMR spectrum. Red spectrum is highly absorbed by vegetation and this is useful to delineate between soil and vegetation.

red1 = extract_band(landsat, [4])
red1

Band ratio technique

A recent remote sensing technique to extract coastline is Band Ratio. In this technique the DN values of bands are divided to create a binary raster. NIR, Red and Blue bands were used for creating the binary raster. NIR band is selected as it is able to delineate water-land boundary, Red band is important for vegetation and water content and Blue band has high reflectance in water bodies. In the current study the following band ratio formula was used:

                                 Blue>NIR & Blue>Red

The output of this formula is a binary raster which has 0 and 1 pixel value. 1 represents water (white pixels) and 0 represents land (black pixels) i.e. bare surface, vegetation etc.

binary1 = arcgis.raster.functions.raster_calculator([nir1, blue1, red1],  
                                                    ['nir', 'blue', 'red'],
                                                    "(blue>nir)&(blue>red)", 
                                                    extent_type='FirstOf', 
                                                    cellsize_type='FirstOf')
binary1

For vizualisation, stretch raster function was used to stretch function the histogram. After applying stretch the pixel values were changed to 0 and 255 where 0 represents land and 255 represents water.

dra1 = arcgis.raster.functions.stretch(binary1, 
                                      stretch_type='MinMax', 
                                      dra=True)
dra1

Postprocessing of results

Clip out extra area

The stretched binary raster was clipped out using aoi geometry to clip out the extra area and to minimize the processing time. Clip raster function was used.

clip_diff1 = clip(dra1, aoi_geom)

water_ras = greater_than([clip_diff1, 0], extent_type='FirstOf', cellsize_type='FirstOf')
water_raster = water_ras.save('water_ras'+str(datetime.now().microsecond), gis=gis2)
water_ras 

dra2 = arcgis.raster.functions.stretch(water_ras, 
                                      stretch_type='MinMax', 
                                      dra=True)
dra2_ras = dra2.save("dra_raster"+str(datetime.now().microsecond), gis=gis2)
dra2

Get the data in feature layer

The clipped binary raster was converted to polygon for coastline extraction. convert_raster_to_feature function was used for conversion.

b2_poly = convert_raster_to_feature(dra2_ras.layers[0], 
                                    field='Value', 
                                    output_type='Polygon', 
                                    simplify=True, 
                                    output_name='coastline_poly'+str(datetime.now().microsecond), 
                                    gis=gis2)

Coastline extraction

Filter out the water polygons

The feature layer was converted to dataframe using query and the polygons were ordered on the basis on Shape_Area column.

## Create dataframe from feature layer and get water polygons
dfm1 = b2_poly.layers[0].query('gridcode=255').sdf 

## Convert dataframe to feature layer
water_poly = gis2.content.import_data(dfm1, title='coast_poly_full'+str(datetime.now().microsecond))
dfm1

Extract the coastline polygon

The polygon with largest area represents the coastline.The polygon with highest value in Shape_Areacolumn was extracted by referring the df. The new dataframe was converted to feature layer using import_data function.

## Get the polygon with largest area as it will represent the coastline
df = water_poly.layers[0].query().sdf
df['MYAREA'] = df.SHAPE.geom.area
dfm5 = df[df['MYAREA']==df['MYAREA'].max()]
coast_poly = gis2.content.import_data(dfm5, title='coast_poly'+str(datetime.now().microsecond))

Convert coastline polygons to line

The coastline polygon was converted to line using PolygonToLine arcpy function.

dc_df = pd.DataFrame.spatial.from_layer(coast_poly.layers[0])
display(dc_df.head())

# Create empty dataframe
df = pd.DataFrame()

# Create polyline
df['SHAPE'] = dc_df.SHAPE.geom.boundary()
coast_polyline = df.spatial.to_featurelayer('coastline_polyline'+str(datetime.now().microsecond))

Remove noise

A mask was created using the clipped raster, it represents the extent of raster. The mask will be used to remove the extra lines.

#Create a mask
mask = greater_than_equal([clip_diff1, -1], extent_type='FirstOf', cellsize_type='FirstOf')
mask_ras = mask.save('mask'+str(datetime.now().microsecond), gis=gis2)

# Create a mask to remove noise which covers whole extent of clipped raster
aoi2 = convert_raster_to_feature(mask_ras.layers[0], 
                                field='Value', 
                                output_type='Polygon', 
                                simplify=True, 
                                output_name='aoi2'+str(datetime.now().microsecond), 
                                gis=gis2)

The coastline layer which we got in the previous step has some extra lines representing the boundary of aoi layer. To remove the extra lines a negative buffer of 10 meter was created using create_buffer function.

### Create neg_buffer to remove the noise
neg_buffer = arcgis.create_buffers(aoi2, 
                                   distances=[-10],
                                   units='Meters', 
                                   output_name='aoi_neg_buf'+str(datetime.now().microsecond), 
                                   gis=gis2)

Extra lines were masked out using the negative buffer boundary.

coastline_final = overlay_layers(coast_polyline.layers[0], neg_buffer.layers[0], output_name='extracted_coastline'+str(datetime.now().microsecond))

Extracted Coastlines

m9= gis2.map('USA')
m9

m9.center = [29.058093925550384, -89.19558942966272]
m9.zoom = 14

m9.content.add(landsat)
m9.content.add(coastline_final)

## Create dataframe for visualization
dfm_1= coastline_final.layers[0].query(out_fields="*").sdf

## add the dataframe to map widget
dfm_1.spatial.plot(map_widget=m9,
                renderer_type='s', 
                pallete='spring',
                alpha=1
                )

True

Coastline for whole USA were extracted with the above workflow using Landsat-8 imagery. To see the results, click here

Conclusion

Band Ratio technique is a efficient method, which gives highly accurate results with less processing time and it is able to cover both temporal and spatial aspects of coastline changes. Band Ratio technique is a easy to calculate method which gives highly accurate results with less processing time. The workflow can be applied on any area using multispectral imagery i.e. Landsat-8, Sentinel-2, etc.

Literature resources