Boat Lights from the VIIRS: Shedding some light on night-time fishing activity in the Indo-Malayan waters

Light emitted at night due to human activities can easily be seen from space and detected by satellite. Many compelling pictures of the Earth at night have been taken revealing cities, roads, and conurbations. Light emitted at sea has, however, received less attention. Applications relating to fisheries and their management might prove particularly useful. NOAA/NASA and partners have recently developed a boat detection product derived from satellite imagery at high temporal and spatial resolution and the data can be downloaded for free from here as either kml/csv files and processed in R. 

For many years (since 1972 to be precise), global low light night imaging was collected by the US Air Force's Defence Meteorological Satellite Programmes Operational Linescan System (DMSP-OLS). While night-time lights collected from DMSP-OLS has major shortcomings, there were no reasonable alternative products for mapping night-lights of lit fishing boats. This all changed in 2011 with the launch of NASA/NOAA Suomi National Polar Partnership satellite equipped with the Visible Infrared Radiometer Suite (VIIRS) as the primary imaging sensor. You can find out many of the advantages the new VIIRS sensor has over the DMPS-OLS here. One of the most important improvements is that the 'Day/Night band' on the VIIRS that is used for night-time light detection has a detection footprint 45 times smaller than that on the DMSP, which ultimately means it can detect light sources a lot more smaller than what the DMPS-OLS could! The VIIRS also has a better dynamic range and lower detection limits which means it can identify and discern more sources of light at night with less saturation when it gets too bright with surrounding lights. These two improvements are key in enabling better detection of boat lights at night.

Commercial fishermen the world over use very bright light to attract squid and small pelagics to their nets, for example, and these are visible from space. Here we use the NOAA boat light data to animate nightly boat lights during February 2015 (cumulative for entire night) for a part of south-east Asia including Indonesia, Malaysia and the northern tip of Australia. We think that many of these lights reflect fishing activity, but more work needs to be done to detect and isolate shipping etc. A concerted fishers' knowledge/consultation process would be essential for enabling us to understand these data properly. We think it would it would be worth the effort, though, and could yield fascinating insights into fleet behavior and the species targeted.

As usual, we provide the R code below on how we produced the above animation of the boat lights detection product.

##### Code to read and produce animation of boat night light detection data from NGDC site #####
 
# Load required libraries
 
library(RCurl)
library(rvest)
library(dplyr)
library(lubridate)
library(zoo)
library(ggplot2)
library(rgdal)
library(grid)
library(animation)
library(gridExtra)
library(extrafont)
 
# Set working directory
 
setwd("D:/Projects/BoatLights")
 
# Load fonts
 
loadfonts()
 
# Set date range for csv files to be downloaded from NGDC site
 
date.range <- seq.Date(from=as.Date('2015-02-01'),as.Date('2015-02-28'),by='1 day')
 
# Get dates into relevant format
 
date <- as.character(date.range)
date1 <- gsub("-","",date)
 
 
# Get urls to scrape boat light csv files off  NGDC site
 
urls <- sprintf((paste("http://mapserver.ngdc.noaa.gov/viirs_data/indo_boats///",c(date1),"/VNB_npp_d",c(date1),
        "_ssec_v12.csv",sep="")))
 
# Read boat light data files in csv format and bind them
 
boats <- do.call(rbind,lapply(urls,read.csv))
 
# Extract useful columns
 
boats.df <- boats[-c(18,19)]
 
# Reformat date and add other date and time columns
 
boats.df$Date <- as.POSIXct(boats.df$Date_Mscan, format="%Y-%m-%d %H:%M:%S",tz="GMT")
boats.df$Dates <- as.Date(boats.df$Date)
boats.df$Year <- as.numeric(as.POSIXlt(boats.df$Date)$year+1900)
boats.df$Month <- as.numeric(as.POSIXlt(boats.df$Date)$mon+1)
boats.df$Monthf <- factor(boats.df$Month,levels=as.character(1:12),
                 labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
boats.df$Weekday = as.numeric(format(as.POSIXlt(boats.df$Date),"%u"))
boats.df$Yearmonth <- as.yearmon(boats.df$Date)
boats.df$Yearmonthf <- factor(boats.df$Yearmonth)
boats.df$Week <- as.numeric(format(as.POSIXlt(boats.df$Date),"%W"))
boats.df$Day <- strftime(boats.df$Date,'%A')
boats.df$Jday <- strptime(boats.df$Date, "%Y-%m-%d")$yday+1
boats.df$Hour <- as.numeric(format(strptime(boats.df$Date, format = "%Y-%m-%d %H:%M"),format = "%H"))
 
# Extract only detections that are classified as 'boat' from quality flag for nightboat detection algorithm
 
boats.df1 <- subset(boats.df,QF_Detect==1)
 
# Import shapefiles
 
world <- readOGR(dsn="D:/GeoSpatial/Data/World_Boundaries", layer="TM_WORLD_BORDERS-0.3") # World boundary shp's download from www.thematicmapping.org
 
eez <- readOGR(dsn="D:/Projects/BoatLights", layer="EEZs") # EEZ boundary shp's downloaded from www.marineregions.org
 
mpa <- readOGR(dsn="D:/Projects/BoatLights", layer="MPAs") # MPA boundary shp's downloaded from CT Atlas website - www.ctlas.reefbase.org
 
# Fortify layers for ggplot mapping
 
world.df <- fortify(world)
 
eez.df <- fortify(eez)
 
mpa.df <- fortify(mpa)
 
# Set theme options for map
 
theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_rect(fill="black"),
                         panel.border = element_rect(colour="black"),
                         legend.text = element_text(size = 8,face="bold",family="Myriad Pro"),
                         legend.position = "right",
                         legend.key.width=unit(0.75,"line"),
                         legend.key.height=unit(0.75,"line"),
                         legend.direction="horizontal",
                         axis.line = element_blank(),
                         axis.text.x = element_blank(),
                         axis.text.y = element_blank(),
                         axis.ticks = element_blank(),
                         axis.title.x = element_blank(),
                         plot.title = element_text(size=20,face="bold",family="Graph Black",hjust=0,vjust=-0.5),
                         axis.title.y = element_blank()))
 
 
# Run loop to produce boat light maps for each night in February 2015
 
for (i in date)
 
{
 
 
b <- ggplot(subset(boats.df1,Dates ==i),aes(Lon_DNB,Lat_DNB))+
     geom_polygon(data=mpa.df, aes(long,lat, group=group,fill="MPA"),colour=NA,alpha=0.5)+
     geom_point(color = "gold", size = 0.5, alpha = 0.5) +
     geom_polygon(data=eez.df, aes(long,lat, group=group,colour="EEZ"),fill=NA,size=0.25,linetype="dashed")+
     geom_polygon(data=world.df, aes(long,lat, group=group),fill="black",colour="darkgray",size=0.35)+
     theme_bw()+ ylab("")+ xlab("")+
     ggtitle((paste("Boat Lights detected from VIIRS on ",i,"\n"))) +
     scale_fill_manual(name="",values=c(MPA="forestgreen"))+
     scale_colour_manual(name="",values=c(EEZ="dimgray"))+
     coord_equal(xlim=c(95,140),ylim=c(-15,9))+
     theme_opts
 
b
 
bg <- arrangeGrob(b, sub = textGrob("Source: VIIRS Boat Detection (VBD) Data (PROTOTYPE) - NGDC (2015)\n", hjust = -0.33,
      vjust=-0.3,gp = gpar(fontface = "bold", fontsize = 8, fontfamily="Tw Cen MT Condensed")))
 
bg
 
ggsave(bg,file=paste("D:/Projects/BoatLights/boatlights_",i,".png",sep=""),width=8,height=5,unit="in",type="cairo-png")
 
}
 
# Set path for Image Magick convert program 
 
path.to.convert <- paste0(shortPathName("C:/Program Files/ImageMagick-6.9.0-Q16)/"),"convert.exe")
ani.options(convert=path.to.convert)
 
 
# Convert png files to GIF animation
 
files = sprintf('boatlights_%s.png', c(date))
im.convert(files, output = 'BoatLights_Feb2015.gif')

Created by Pretty R at inside-R.org

Blowing in the wind - How many times must we keep experiencing the haze before the Governments do something about it?

Key points of post

• Visibility in Penang is affected by the seasonal monsoon winds
• Visibility has got worse since 2001.
• The governments of south-east Asia must urgently solve this potentially serious health issue.

Every year South-East Asians suffer from the effects of fires burning in the forests of Indonesia. The more intense fires are lit deliberately to clear land for palm oil plantations. It is illegal to use fire to clear land in Indonesia, but the practice is unfortunately widespread, leading to habitat destruction and loss of biodiversity. A satellite image captured by the Moderate Resolution Imaging Spectroradiomater (MODIS) from NASA below for March 2014 illustrates when the fires were particularly bad (fires shown in red). 

Source: NASA Earth Observatory (2014) - http://1.usa.gov/1kqfdYA

Most fires occur during the dry season between April and October. The fires in 2014 started in February, and actually forced Indonesia to declare a state of emergency because of poor air quality. 

But what about the situation for us in Penang ? Our island is subject to two seasonal monsoon winds: the north-east and the south-west. These patterns are clear from the wind-rose diagram below which shows the average monthly wind-speeds and directions for the period 1949-2014. 

Given the location of Indonesia relative to Penang, it follows that visibility will be worst during the south-west monsoon which blows hardest, and most consistently, between May and August each year. Conversely we might expect visibility to be best during the north-east monsoon. This is indeed exactly what we observe. Average visibility readings from Bayan Lepas International Airport by week and year 2001-2014 are plotted below. The white squares represent good visibility and cleaner air; the dark squares poor visibility and high levels of air pollution.

It is clear that visibility has gotten worse since 2001, ie. the numbers of white squares have decreased. The pattern is also strongly seasonal which reflects: (a) the monsoons; and (b) the fact that fires are used more in Indonesia’s dry season. The relatively long periods of good visibility (>10 km's) that Penangites must have enjoyed especially during January, November and December seem now to be a thing of the past. I’m (Doug) a keen hiker regularly gazing down on the Island from the heights of Tea Station 5, 89, and Penang Hill. I’ve only been here since January 2012 and, judging from these data, I guess I've probably never seen a really good view!

Let’s hope the Government of Malaysia / Penang will continue to apply pressure and work with the Indonesians and other countries in the region to solve this terrible problem.

If you're keen on producing similar plots like the ones we've included in this blog post, you can use the code provided below. The wind-rose plot was produced using the 'openair' package for R which is an open-source air pollution analysis tool. There are a number of useful plots that can be produced using this package to show magnitude and direction of wind / air pollutants which you can further explore on your own. We also provide you with the code to produce the 'heatmap' plot of the visibility data which was pulled into R from Weather Underground (using the 'weatherData' package).  

# Load libraries
 
library(ggplot2)
library(RColorBrewer)
library(openair)
library(weatherData)
library(mgcv)
library(scales)
library(plyr)
library(reshape2)
library(circular)
library(gridExtra)
library(lubridate)
library(weathermetrics)
library(zoo)
 
# Setting work directory
 
setwd("d:\\ClimData")
 
# Reading and reformatting raw daily data downloaded from NCDC
 
dat<-read.table("2422706962434dat.txt",header=TRUE,fill=TRUE,na.strings=c("*","**","***","****","*****","******","0.00T*****"))
 
colnames(dat)<-tolower(colnames(dat))
 
Sys.setenv(TZ = "UTC")
dat$date <- as.POSIXct(strptime(dat$yr..modahrmn,format="%Y%m%d%H%M"))  + 8 * 60 * 60
 
dat$dir[dat$dir == 990.0] <- NA 
 
# Convert windspeed in miles/hour to metres/second
 
dat$wspd <- (dat$spd)*0.44704
 
# Create columns for year, month and dat
 
dat$year <- as.numeric(format(dat$date,"%Y"))
 
dat$month <- as.numeric(format(dat$date,"%m"))
 
dat$day <- as.numeric(format(dat$date,"%d"))
 
# Subset year of interest
 
datsub <- subset(dat, year >= 1949 )
 
# Remove NA's
 
datsub$dir[is.na(datsub$dir)] <- 0
 
# Plot and produce png file of Wind Rose
 
png(filename = "Monthly_WindRose_Penang_1949-2014.png",height=10,width=10,
    bg = "white",units='in', res = 600, family = "", restoreConsole = TRUE,
    type = "cairo-png")
 
windRose(datsub,ws="wspd",wd="dir",bias.corr = TRUE,border=TRUE,type = "month",width = 0.5,grid.line = 10,
         statistic = "prop.mean",offset = 10, paddle =FALSE,cols="hue",annotate=FALSE,auto.text=FALSE,
         ws.int=2,breaks=c(0,2,4,6,8,10,12),key = TRUE, key.footer = "(meter/second)", key.position = "bottom",na.action=NULL,
         key.header = "Wind Speed",main="       Average Monthly Wind Rose plots for Penang (Bayan Station)\n from years 1949 - 2014 \n",font.main=1,
         sub="\nData source: Integrated Surface Database (ISD) - National Climatic Data Centre (NCDC)",font.sub=2)
 
dev.off()
 
###########################################
#### Produce visibility 'heatmap' plot ####
###########################################
 
### Find required 4 letter station code using function below
 
# Example getStationCode("Penang") or getStationCode("Honiara")
 
# Information can also be acquired from the following link - http://weather.rap.ucar.edu/surface/stations.txt
 
getStationCode("Butterworth")
 
# Use station code below to get required plot data and parameters
 
station.id="WMKP"
 
s<-getStationCode(station.id)
 
s1<-(strsplit(s,split= " "))[[1]]
 
station.name<-paste(s1[4],s1[5])
 
### Getting summarized weather data for WU
 
ws2001<-getSummarizedWeather(station.id, "2001-01-01", "2001-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2002<-getSummarizedWeather(station.id, "2002-01-01", "2002-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2003<-getSummarizedWeather(station.id, "2003-01-01", "2003-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2004<-getSummarizedWeather(station.id, "2004-01-01", "2004-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2005<-getSummarizedWeather(station.id, "2005-01-01", "2005-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2006<-getSummarizedWeather(station.id, "2006-01-01", "2006-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2007<-getSummarizedWeather(station.id, "2007-01-01", "2007-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2008<-getSummarizedWeather(station.id, "2008-01-01", "2008-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2009<-getSummarizedWeather(station.id, "2009-01-01", "2009-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2010<-getSummarizedWeather(station.id, "2010-01-01", "2010-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2011<-getSummarizedWeather(station.id, "2011-01-01", "2011-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2012<-getSummarizedWeather(station.id, "2012-01-01", "2012-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2013<-getSummarizedWeather(station.id, "2013-01-01", "2013-12-31", opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
ws2014<-getSummarizedWeather(station.id, "2014-01-01", Sys.Date(), opt_custom_columns=T,

custom_columns=c(1:23), opt_verbose=T)
 
ws<-rbind(ws2001,ws2002,ws2003,ws2004,ws2005,ws2006,ws2007,ws2008,ws2009,ws2010,ws2011,ws2012,ws2013,ws2014) 
 
summary(ws)
 
str(ws)
 
### List of variables from Wunderground ###
 
# [1] "MYT"                        "Max_TemperatureC"           "Mean_TemperatureC"          "Min_TemperatureC"          
# [5] "Dew_PointC"                 "MeanDew_PointC"             "Min_DewpointC"              "Max_Humidity"              
# [9] "Mean_Humidity"              "Min_Humidity"               "Max_Sea_Level_PressurehPa"  "Mean_Sea_Level_PressurehPa"
# [13] "Min_Sea_Level_PressurehPa"  "Max_VisibilityKm"           "Mean_VisibilityKm"          "Min_VisibilitykM"          
# [17] "Max_Wind_SpeedKm_h"         "Mean_Wind_SpeedKm_h"        "Max_Gust_SpeedKm_h"         "Precipitationmm"           
# [21] "CloudCover"                 "Events"                     "WindDirDegrees"            
 
colnames(ws)<-c("date","date1","maxtemp","meantemp","mintemp","dewp","meandewp","maxdewp","maxhum","meanhum","minhum","maxslp","meanslp",
                "minslp","maxvsb","meanvsb","minvsb","maxwspd","meanwspd","maxgust","prcp","cc","events","wd")
 
## Adding date columns

ws$dates <- as.Date(ws$date)
ws$year <- as.numeric(as.POSIXlt(ws$dates)$year+1900)
ws$month <- as.numeric(as.POSIXlt(ws$dates)$mon+1)
ws$monthf <- factor(ws$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)
ws$weekday <- as.POSIXlt(ws$dates)$wday
ws$weekdayf <- factor(ws$weekday,levels=rev(0:6),labels=rev(c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")),ordered=TRUE)
ws$yearmonth <- as.yearmon(ws$date)
ws$yearmonthf <- factor(ws$yearmonth)
ws$week <- as.numeric(format(as.Date(ws$dates),"%W"))
ws$weekf<- factor(ws$week)
ws$jday<-yday(ws$dates)
 
# Define colour palette
 
col<-c("black","grey10","grey20","grey50","white")
 
# Plot 'heatmap'
 
v <- ggplot(data=ws,aes(x=week,y=year,fill=meanvsb))+ geom_tile(colour="black",size=0.65)+ theme_bw()+
     scale_fill_gradientn(colours=col,name="Visibility\n(km)\n")+coord_equal(ratio=1)+
     ylab("YEAR\n")+
     xlab("\nWEEK OF YEAR\n\nSource: Weather Underground (2014)")+
     scale_y_continuous(expand = c(0,0),breaks = seq(2000, 2015, 1)) + scale_x_discrete(expand = c(0,0),breaks = seq(0,52,2))+
     ggtitle("Average weekly visibility in Penang\n")+
     theme(panel.background=element_rect(fill="transparent"),
        panel.border=element_blank(),
        axis.title.y=element_text(size=10,colour="grey20"),
        axis.title.x=element_text(size=10,colour="grey20"),
        axis.text.y=element_text(size=10,colour="grey20",face="bold"),
        axis.text.x=element_text(size=10,colour="grey20",face="bold"),
        plot.title = element_text(lineheight=1.2, face="bold",size = 14, colour = "grey20"),
        panel.grid.major = element_blank(),
        legend.key.width=unit(c(0.2,0.2),"in"))
 
v
 
# Save plot to png
 
ggsave(v, file="Penang_Weekly_Average_Visibility.png", width=15, height=5,dpi=400,type = "cairo-png")

Created by Pretty R at inside-R.org