Loading ship detections via Sarscape processes

wd <- "C:/Users/Robert/Desktop/shipScore/analytics/shiny_ship_detection/import/S1A_IW_GRDH_1SDV_20191223T142452_20191223T142517_030477_037D47_221F.SAFE/"

darkobj_pt <- shapefile(paste0(wd,"sentinel1_SENTINEL-1_130_20191223_142452805_IW_SIW_A_VV_gr_geo_ship_mask.shp"))
darkobj_poly <- shapefile(paste0(wd,"ship_vectors_multi.shp"))

Loading AIS signals via cassandra output

# AIS data load
ais_signals_uncorrected <- read.csv('C:\\Users\\Robert\\Desktop\\ISR-Maritime-Analytics\\analytics\\tests\\data-1579797017584.csv', header = TRUE, stringsAsFactors=F)
coordinates(ais_signals_uncorrected) <- ~lng+lat

Preparing test shorelines using GIS

• Shoreline data retrieved from https://dnc.nga.mil/PrototypeGSD.php
  • Retrieval date unknown, circa 2004.
  • Data is said to be high water line.
• Inside ArcGIS
  • Used Feature to Polygon with the footprint/envelope and shoreline polyline as inputs.
  • Manually removed ocean features inside an Edit.
  • Perform a ‘Select By Location’ and delete the inverse (using switch selection inside attribute table).
    • Target layer: shoreline polygon output
    • Source layer: footprint/envelope
    • Spatial selection method: ’are within (Clementini) the source layer feature.
  • Visual spot check against 1/9/2019 Digital Globe Satellite Imagery
    • Various infill and “made land” not included from between approx 2004 and 2019.
    • In the United Arab Emirites, there are numerous such examples of recent development.

Loading test shorelines

• Retrieved from http://www.soest.hawaii.edu/pwessel/gshhg/

wd <- "C:/Users/Robert/Desktop/ISR-Maritime-Analytics/analytics/shoreline_data/"
shoreline_poly <- shapefile(paste0(wd,"nga_zone10_example.shp"))

Apply secondary shoreline filter

# NAs are points not on land
darkobj_pt_filtered<-darkobj_pt[!complete.cases(over(darkobj_pt, shoreline_poly)),]

Data Visualization

leaflet(ais_signals_uncorrected) %>% 
  addProviderTiles("CartoDB.Positron") %>%
  setView(55,26, zoom = 10) %>%
  # Currently representing AIS signals
  addCircleMarkers(label=~mmsi,
                   weight = 3, 
                   radius=10, 
                   color="#ffa500")  %>%
  # Currently representing dark object points
  addMarkers(data = darkobj_pt_filtered,
             icon = list(
            iconUrl = 'http://pngimg.com/uploads/target/target_PNG38.png',
            iconSize = c(25, 25)
           ))  %>%
  # Currently representing dark object polygon shapes
  addPolygons(data=darkobj_poly,
              col = 'red',
              stroke = TRUE, 
              #fillOpacity = 0.3, 
              smoothFactor = 0.5) %>%
  addPolygons(data=shoreline_poly,
              col = 'green',
              stroke = FALSE, 
              #fillOpacity = 0.3, 
              smoothFactor = 0.5)

Discussion

• We notice an improvement in the filtering of false positives from the low resolution shoreline output from sarscape. This was expected and pending an improved shoreline feature extraction system, ought to be applied to sarscape’s methods prior to this point.
• The shoreline has substantially shifted, leaving many false positives along made land. Further details here: https://en.wikipedia.org/wiki/Land_reclamation_in_the_United_Arab_Emirates.
• AIS signals themselves appear within the land mask.
• An improved product may be available for purchase at approx $200 from https://shoreline.noaa.gov/data/datasheets/wvs.html.

Data tidying

` - Remove “distant” points greater than 30 seconds.

ais_signals_uncorrected_30s <- subset(ais_signals_uncorrected, delta_secs < 30)

#http://www.sthda.com/english/articles/32-r-graphics-essentials/133-plot-one-variable-frequency-graph-density-distribution-and-more/
     
# Change outline and fill colors by groups ("sex")
# Use a custom palette
ggdensity(ais_signals_uncorrected_30s@data, x = "delta_secs",
          fill = "#0073C2FF", color = "#0073C2FF",
          add = "mean", rug = TRUE)

• Conforming projections in UTM for geoprocessing accuracy

# calculate UTM Zone
long2UTM <- function(long) {
    (floor((long + 180)/6) %% 60) + 1
}
centroid <- shoreline_poly@polygons[[1]]@labpt
if (centroid[2] > 0){ hemisphere = 'N'}else{hemisphere = 'S'}

utm_zone <- paste0("+proj=utm +zone=",paste0(long2UTM(centroid[1]) , hemisphere)," +datum=WGS84 +units=km")

proj4string(ais_signals_uncorrected_30s) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" )

ais_signals_uncorrected_30s_UTM <- spTransform(ais_signals_uncorrected_30s, CRS(utm_zone))
shoreline_poly_UTM <- spTransform(shoreline_poly, CRS(utm_zone))

g = rep(NA, dim(ais_signals_uncorrected_30s_UTM)[1])
for(i in 1:dim(ais_signals_uncorrected_30s_UTM)[1]){
  g[i] = gDistance(ais_signals_uncorrected_30s_UTM[i,],shoreline_poly_UTM)
}
ais_signals_uncorrected_30s_UTM@data$dShore <- g

ggdensity(ais_signals_uncorrected_30s_UTM@data, x = "dShore",
          fill = "#0073C2FF", color = "#0073C2FF",
          add = "mean", rug = TRUE)

• Remove signals within 3 km of shoreline and reproject UTM dataset into decimal degrees for mapping.

ais_signals_uncorrected_30s_UTM_3km <- subset(ais_signals_uncorrected_30s_UTM, dShore > 3)
ais_signals_uncorrected_30s_wgs84_3km <- ais_signals_uncorrected_30s_UTM_3km
ais_signals_uncorrected_30s_wgs84_3km <- spTransform(ais_signals_uncorrected_30s_wgs84_3km, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

This leaves 11 candidate signals out of an original 326 vessels`

Visualization

leaflet(ais_signals_uncorrected_30s_wgs84_3km) %>% 
  addProviderTiles("CartoDB.Positron") %>%
  setView(55,26, zoom = 8) %>%
  # Currently representing AIS signals
  addCircleMarkers(label=~mmsi,
                   weight = 3, 
                   radius=20, 
                   color="#ffa500")  %>%
  # Currently representing dark object points
  addMarkers(data = darkobj_pt_filtered,
             icon = list(
            iconUrl = 'http://pngimg.com/uploads/target/target_PNG38.png',
            iconSize = c(25, 25)
           ))  %>%
  # Currently representing dark object polygon shapes
  addPolygons(data=darkobj_poly,
              col = 'red',
              stroke = TRUE, 
              #fillOpacity = 0.3, 
              smoothFactor = 0.5) %>%
  addPolygons(data=shoreline_poly,
              col = 'green',
              stroke = FALSE, 
              #fillOpacity = 0.3, 
              smoothFactor = 0.5)