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

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.

image.png

Neccessary Imports

Input
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_feature_to_raster, convert_raster_to_feature
from arcgis.geoanalytics.manage_data import clip_layer
from arcgis.raster.functions import equal_to, greater_than, clip, apply, extract_band, stretch

Connect to your GIS

Input
from arcgis import GIS
gis =  GIS('home')
gis2 = GIS('https://pythonapi.playground.esri.com/portal', 'arcgis_python', 'amazing_arcgis_123')

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.

Input
landsat_item = gis.content.search('title:Multispectral Landsat owner:esri_livingatlas tags:Multitemporal, imagery, landsat 8, temporal, MS', 'Imagery Layer', outside_org=True)[0]
landsat = landsat_item.layers[0]
landsat_item
Output
Multispectral Landsat
Landsat multispectral and multitemporal imagery with on-the-fly renderings and indices for visualization and analysis. The Landsat 8 imagery in this layer is updated daily and is directly sourced from the Landsat on AWS collection.Imagery Layer by esri_livingatlas
Last Modified: May 14, 2020
0 comments, 1 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.

Input
aoi = gis.content.search('usa_coast_buff_f', 'feature layer')[0]
aoi
Output
usa_coast_buff_f
usa_coast_buffFeature Layer Collection by api_data_owner
Last Modified: September 28, 2020
0 comments, 9 views

Prepare data for analysis

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

Input
m = gis2.map('USA', 4)
m

Input
m.add_layer(aoi)
m.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
 'xmin': -10313326.840304691,
 'ymin': 3355270.9128910713,
 'xmax': -9711614.553643815,
 'ymax': 3599869.403403623}

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.

Input
dfm = aoi.layers[0].query(out_fields="OBJECTID, type, code").sdf
dfm
Output
OBJECTID type code SHAPE
0 1 islands 1 {"rings": [[[-17328342.0256, 2126164.6003], [-...
1 2 eastern_coast 1 {"rings": [[[-10796721.2644, 2984284.5892], [-...
2 3 alaska 1 {"rings": [[[-17408919.7589, 11517244.7315], [...
3 4 western_coast 1 {"rings": [[[-13019578.4789, 3857089.5857], [-...

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.

Input
aoi_layer = aoi.layers[0]
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.

Input
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
Output
OBJECTID AcquisitionDate GroupName CloudCover DayOfYear Shape_Length Shape_Area SHAPE
0 2394944 2018-12-21 16:51:19 LC80250412018355LGN00_MTL 0.0000 355 8.406314e+05 4.412466e+10 {"rings": [[[-10505837.9096, 3260395.169700000...
1 2973331 2017-03-04 16:07:36 LC80180402017063LGN00_MTL 0.0000 63 8.520258e+05 4.532983e+10 {"rings": [[[-9261778.5607, 3442875.9834999964...
2 3146484 2019-08-30 15:34:41 LC80130342019242LGN00_MTL 0.0000 242 9.400965e+05 5.519010e+10 {"rings": [[[-8120177.4793, 4595814.804499999]...
3 1972432 2016-10-24 15:34:50 LC80130342016298LGN00_MTL 0.0000 298 9.411144e+05 5.530684e+10 {"rings": [[[-8122927.2599, 4596153.234999999]...
4 1985184 2017-10-17 16:38:35 LC80230392017290LGN00_MTL 0.0000 290 8.641411e+05 4.664550e+10 {"rings": [[[-10094115.3386, 3565121.854099996...
... ... ... ... ... ... ... ... ...
782 2769986 2016-05-07 16:38:16 LC80230402016128LGN01_MTL 0.0293 128 8.519905e+05 4.532605e+10 {"rings": [[[-10123608.5317, 3442801.410499997...
783 2840922 2017-04-09 15:40:58 LC80140362017099LGN00_MTL 0.0296 99 9.065239e+05 5.131708e+10 {"rings": [[[-8394085.965, 4199193.041100003],...
784 2758546 2016-04-11 16:01:36 LC80170412016102LGN01_MTL 0.0297 102 8.407949e+05 4.414217e+10 {"rings": [[[-9132542.3423, 3260274.1598000005...
785 2747763 2016-03-06 16:26:12 LC80210402016066LGN01_MTL 0.0299 66 8.520226e+05 4.533050e+10 {"rings": [[[-9777993.9307, 3442750.3139000013...
786 1971290 2017-06-14 15:26:01 LC80120292017165LGN00_MTL 0.0299 165 1.047326e+06 6.850239e+10 {"rings": [[[-7672954.216399999, 5657896.10689...

787 rows × 8 columns

Input
## 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.

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

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.

Input
nir1 = extract_band(landsat, [5])
nir1
Output

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

Input
blue1 = extract_band(landsat, [2])
blue1
Output

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.

Input
red1 = extract_band(landsat, [4])
red1
Output

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.

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

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.

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

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.

Input
clip_diff1 = clip(dra1, aoi_geom)
clip_diff_ras1 = clip_diff1.save("b_poly1oct", gis=gis2)
Input
water_ras = greater_than([clip_diff1, 0], extent_type='FirstOf', cellsize_type='FirstOf')
water_ras 

dra2 = arcgis.raster.functions.stretch(water_ras, 
                                      stretch_type='MinMax', 
                                      dra=True)
dra2
Output

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.

Input
b2_poly = convert_raster_to_feature(dra2, 
                                    field='Value', 
                                    output_type='Polygon', 
                                    simplify=True, 
                                    output_name=None, 
                                    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.

Input
## Create dataframe from feature layer and get water polygons
dfm1 = b2_poly.layers[0].query('gridcode=1').sdf 
## Convert dataframe to feature layer
water_poly = gis2.content.import_data(dfm1, title='test_poly')
dfm1
Output
FID Id gridcode Shape__Area Shape__Length SHAPE
0 1 1 1 5.763496e+02 1.124948e+02 {"rings": [[[-10148539.824, 3424382.5147], [-1...
1 2 2 1 5.763496e+02 1.124948e+02 {"rings": [[[-10141969.8248, 3424382.5147], [-...
2 3 3 1 5.763496e+02 1.124948e+02 {"rings": [[[-10141069.8249, 3424382.5147], [-...
3 4 4 1 9.000000e+02 1.200000e+02 {"rings": [[[-10139629.8251, 3424352.5147], [-...
4 7 7 1 5.763496e+02 1.124948e+02 {"rings": [[[-10135069.8257, 3424382.5147], [-...
... ... ... ... ... ... ...
5361 7457 7457 1 1.214914e+03 1.702172e+02 {"rings": [[[-10119349.8277, 3394052.5185], [-...
5362 7480 7480 1 9.000039e+02 1.200000e+02 {"rings": [[[-10116349.8281, 3393812.5186], [-...
5363 7485 7485 1 1.045953e+03 1.726506e+02 {"rings": [[[-10114789.8283, 3393812.5186], [-...
5364 7489 7489 1 1.800000e+03 1.800000e+02 {"rings": [[[-10113439.8284, 3393812.5186], [-...
5365 7491 7491 1 1.444475e+09 2.444402e+06 {"rings": [[[-10113379.8284, 3393812.5186], [-...

5366 rows × 6 columns

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.

Input
## Get the polygon with largest area as it will represent the coastline
df=water_poly.layers[0].query().sdf
dfm5 = df[df['SHAPE__Area']==df['SHAPE__Area'].max()]
coast_poly=gis2.content.import_data(dfm5, title='coastpoly1oct_1')
dfm5
Output
objectid fid id gridcode SHAPE__Length SHAPE__Area SHAPE
5365 5366 7491 7491 1 2.444402e+06 1.444475e+09 {"rings": [[[-10187209.819100002, 3420632.5151...

Convert coastline polygons to line

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

Input
## To get the url which will be used as the input for PolygonToLine function
coast_poly.url
Output
'https://deldevd040.esri.com/server/rest/services/Hosted/aaeb81/FeatureServer'
Input
line = arcpy.management.PolygonToLine("https://deldevd040.esri.com/server/rest/services/Hosted/aaeb81/FeatureServer/0", 
                                      r"C:\data\coastline1.shp",
                                      "IGNORE_NEIGHBORS")
Input
#define location of shapefiles and destination of zipped shapefiles
source = r"C:\data"
dest = r"C:\data\zip_shp"
filename = 'coastline1'
 
#change the current directory
os.chdir(source)
 
shapefile = f"{filename}.shp"
print(shapefile)
 
# create destination directory if it does not exist
if not os.path.exists(dest):
    os.makedirs(dest)

#creates the name for each zipefile based on shapefile root names
dest_file = os.path.join(dest, filename+".zip")
with ZipFile(dest_file, "w") as z: 
    for file in glob.glob(f'{filename}*'):
        z.write(file)
Input
## The output of PolygonToLine function was published for further analysis
data = r"C:\Users\shi10484\Downloads\zip_shp\coastline1.zip"
shpfile = gis2.content.add({}, data)
coastline = shpfile.publish()

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.

Input
# Create a mask to remove noise which covers whole extent of clipped raster
aoi2 = clip_diff1>=0
aoi2

### Create neg_buffer to remove the noise
neg_buffer = arcgis.create_buffers(aoi2, 
                                   distances=[-10],
                                   units='Meters', 
                                   output_name='aoi_neg_buffer10m', 
                                   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.

Input
### Create neg_buffer geom to mask out the noise

aoi1_layer = neg_buffer.layers[0]
aoi1_feature = aoi1_layer.query()
aoi1_geom = aoi1_feature.features[0].geometry
aoi1_geom['spatialReference'] = {'wkid':3857}

Extra lines were masked out using the negative buffer boundary.

Input
coastline_final = arcgis.geoanalytics.manage_data.clip_layer(coastline, 
                                                             neg_buffer, 
                                                             output_name='coastline_f', 
                                                             gis=gis2)

Extracted Coastlines

Input
m9= gis.map('USA', 4)
m9

Input
m9.add_layer(landsat)
Input
## add the dataframe to map widget
dfm5.spatial.plot(map_widget=m9,
                renderer_type='s', 
                pallete=[255, 255, 0, 0],
                alpha=1
                )
Output
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

Literature Source Author
Research Paper Shoreline Change Mapping Using Remote Sensing and GIS Ali Kourosh Niya, Ali Asghar Alesheikh, Mohsen Soltanpor, Mir Masoud Kheirkhahzarkesh
Research Paper Shoreline change assessment using remote sensing and GIS techniques: a case study of the Medjerda delta coast, Tunisia Mourad Louati & Hanen Saidi & Fouad Zargouni

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.