Mapfishing

Configuramos el directorio madre.

setwd("C:/Users/ADMIN/OneDrive/mantas/análisis embarcaciones")

Importing data

library(tidyverse)
library(data.table)



# List all the CSV files in your folder
csv_files <- list.files(path = "boats/fishing", pattern = "*.csv", full.names = TRUE)
csv_files2 <- list.files(path = "boats/water_taxi", pattern = "*.csv", full.names = TRUE)

# Import and merge all the CSV files into one data frame
fishing <- lapply(csv_files, read_csv) %>% rbindlist()
water_taxi <- lapply(csv_files2, read_csv) %>% rbindlist()
#Add column names
colnames(fishing) <- c("id","lat","lon","date","speed")
colnames(water_taxi) <- c("id","lat","lon","date","speed")

Setting up time variables

Using the “Lubridate” library we create two date variables and correct the time zone.

library(lubridate)

fishing$dmy_hm<-dmy_hm(fishing$date,truncated = 1,tz = "UTC")
fishing$dmy_hm<- with_tz(fishing$dmy_hm, tzone = "America/Mexico_City")
fishing$date<-substr(fishing$dmy_hm,1,10)


water_taxi$dmy_hm<-dmy_hm(water_taxi$date,truncated = 1,tz = "UTC")
water_taxi$dmy_hm<- with_tz(water_taxi$dmy_hm, tzone = "America/Mexico_City")
water_taxi$date<-substr(water_taxi$dmy_hm,1,10)

Parsing the dataset with the boat´s metadata

In this step, valuable information like boat´s name, home port and fishing gear is included in the dataset

metadata<-read.csv("metadata.csv",header = TRUE)

metadata<-metadata[,c(2,3,4,6,7,8)]




left_join(fishing,metadata,by="id")->fishing

left_join(water_taxi,metadata,by="id")->water_taxi


fishing$id<-as.factor(fishing$id)

Filter out observations in fishing towns

We need to exclude data recorded inside fishing towns and touristic ports because it will bias the analysis of the areas were mantas and boats interact the most.

## Warning: package 'sf' was built under R version 4.3.3

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

## Reading layer `Boca.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Boca.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.3245 ymin: 20.50978 xmax: -105.3136 ymax: 20.51819
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Chimo.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Chimo.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.5991 ymin: 20.47313 xmax: -105.5869 ymax: 20.47782
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Corrales.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Corrales.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.681 ymin: 20.40531 xmax: -105.6689 ymax: 20.41451
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Lacruz.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Lacruz.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.3827 ymin: 20.74533 xmax: -105.3706 ymax: 20.7523
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Los muertos.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Losmuertos.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.3524 ymin: 20.50608 xmax: -105.3375 ymax: 20.51182
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Marina Vallarta.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Marinavallarta.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.259 ymin: 20.64831 xmax: -105.2355 ymax: 20.66785
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Muelle_muertos.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Muellemuertos.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.2423 ymin: 20.59707 xmax: -105.2387 ymax: 20.60231
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Pizota.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Pizota.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.4935 ymin: 20.48895 xmax: -105.4844 ymax: 20.49411
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Punta de Mita.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Puntamita.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.5206 ymin: 20.76769 xmax: -105.5135 ymax: 20.77293
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Quimixto.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Quimixto.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.3739 ymin: 20.50036 xmax: -105.3619 ymax: 20.50769
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Yelapa.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Yelapa.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.4548 ymin: 20.48829 xmax: -105.4414 ymax: 20.49773
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Mismaloya.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Mismaloya.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.3001 ymin: 20.53034 xmax: -105.2872 ymax: 20.54026
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## Reading layer `Colomitos.kml' from data source 
##   `C:\Users\ADMIN\OneDrive\mantas\análisis embarcaciones\kml\Colomitos.kml' 
##   using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XYZ
## Bounding box:  xmin: -105.3277 ymin: 20.51021 xmax: -105.3236 ymax: 20.51353
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84

## st_as_s2(): dropping Z and/or M coordinate

## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate

#exploring water taxi activity

#exploring fishing activity

#Rutas del water taxi

library(ggmap)

## Warning: package 'ggmap' was built under R version 4.3.3

## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
##   Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
##   OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.

mykey<-"AIzaSyCsouMjaesCVvfrydxBQAT--ysO-xtp0g4"
register_google(key=mykey) 

water_taxi$IDdate<-paste(water_taxi$ID,water_taxi$dmy_hm)


ggmap(get_googlemap(center = 
                      c(lon=-105.4,lat=20.52), 
                    zoom = 11,maptype="satellite"))+
  theme_void()->mapBoca

## ℹ <https://maps.googleapis.com/maps/api/staticmap?center=20.52,-105.4&zoom=11&size=640x640&scale=2&maptype=satellite&key=xxx--ysO-xtp0g4>

mapBoca+
  geom_hex(data=water_taxi,aes(lon,lat),bins=40,alpha=0.5)+
  scale_fill_distiller(palette="Spectral")+
  labs(fill="Detections",title = "Water taxi areas",x="Longitude",y="Latitude")+
  theme_minimal(base_size = 20)->watertaxi_activity


mapBoca+
  geom_point(data=water_taxi%>%
               filter(speed>5),aes(lon,lat,col=speed),alpha=0.8,cex=0.6)+
  scale_color_gradientn(colors = c("white", "yellow", "orange", "red"),
                        values = c(0.0,0.25, 0.5, 0.75, 1))+
  labs(col="Speed (Km/h)",title = "Water taxi routes",x="Longitude",y="Latitude")+
  theme_minimal(base_size = 20)->watertaxi_activity2


mapBoca+
  geom_point(data=water_taxi%>%
               filter(speed>40),aes(lon,lat,col=speed),alpha=0.8,cex=0.5)+
  scale_color_gradientn(colors = c("orange", "red"),
                        values = c(0.0,0.25, 0.5, 0.75, 1))+
  labs(col="Speed (Km/h)",title = "Water taxi routes",x="Longitude",y="Latitude")+
  theme_bw(base_size = 20)->watertaxi_activity3

mapBoca+
  stat_summary_hex(data=
                     water_taxi%>%
               filter(speed>40),
                   aes(lon,lat,z=speed),
                   fun = function(x) mean(x), 
                   bins = 100) +
  scale_fill_gradientn(colors = c("white", "yellow", "orange", "red"),
                        values = c(0.0,0.25, 0.5, 0.75, 1))+
  labs(fill="Speed (Km/h)",y="Latitude",x="Longitude",subtitle = "Water taxi High speed areas (above 40 km/h)")+
  theme_bw(base_size = 18)->watertaxi_activity4


watertaxi_activity

## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_binhex()`).

watertaxi_activity2

watertaxi_activity3

watertaxi_activity4

#ggsave("graphics/231217_water taxi routes in the southern coast_2.png",watertaxi_activity2,dpi=500,height = 12,width = 12)
#ggsave("graphics/231217_water taxi routes in the southern coast_3.png",watertaxi_activity3,dpi=500,height = 12,width = 12)
#ggsave("graphics/231217_water taxi routes in the southern coast_4.png",watertaxi_activity4,dpi=500,height = 12,width = 12)

#Rutas de pesca

#Calculating stopped times of fishing vessels

#######programming a stopped TRUE/FALSE variable###########

The idea is first identify the moments when the boats stop in a regular day. I grouped the observations by the vessel id and date to calculate a categorical variable with two factors called stopped based on the speed pattern of each vessel.

Things considered to calculate this where:

a series of observations with 0 speed values means a moment when the vessel stopped.
2 is the minimum number of consecutive observations with 0 speed values needed to classify them as stopped.
Possible drift and unsignificant movements among a group of stopped observations are considered as stopped moment if a single observation with speed greater than 0 is found among a group of observations with 0 speed value)

fishing %>%
  group_by(id, date) %>% #fishing data is grouped by name of the boat and the date (YMD).
  mutate(stopped = ifelse(speed == 0 & (row_number() <= 2 | lag(speed, 2) >= 0), TRUE, FALSE)) %>% #1
  mutate(stopped = ifelse(stopped == TRUE | (row_number() > 1 & speed > 0 & lag(stopped) == TRUE & lead(speed)  == 0), TRUE, stopped)) %>% #2
  mutate(stopped = ifelse(speed == 0 & lead(stopped) == TRUE, TRUE, stopped)) %>%
  mutate(stopped = ifelse((row_number() == n()) & speed == 0, TRUE, stopped)) %>%
  mutate(stopped = ifelse(speed == 0 & lag(speed) != 0 & lead(speed) != 0, FALSE, stopped)) %>%
  mutate(stopped = ifelse(speed == 0 & lag(stopped) == TRUE, TRUE, stopped)) %>%
  mutate(stopped = ifelse(speed > 0 & lag(speed) == 0 & lead(speed) == 0, TRUE, stopped)) %>%
  mutate(stopped = ifelse(speed == 0 & (lag(stopped) == TRUE | lead(stopped) == TRUE), TRUE, stopped)) %>%
  mutate(stopped = ifelse(stopped == FALSE & lag(stopped) == TRUE & lead(stopped) == TRUE, TRUE, stopped))%>%
  mutate(stopped = ifelse(is.na(stopped), FALSE, stopped))%>%
  mutate(stopped = ifelse(lead(dmy_hm) - dmy_hm >= hours(2), FALSE, stopped))%>%
  mutate(stopped = ifelse(is.na(stopped), FALSE, stopped))->fishing

fishing data is grouped by name of the boat and the date (YMD).

(1 Creamos la variable stopped indicando que si la velocidad es de 0 y el número de hilera es menor o = a 2 y la velocidad anterior es mayor o igual a 0. Esto captura los casos cuando un bote se encuentra detenido al inicio de cada día y cuando se detiene despues de un periodo de movimiento.

The boat is stopped if the previous condition was true or if the row number is greater than one and the speed is greater than zero and the previous value of stopped was true and the next value of speed is zero. This condition captures the cases where the boat moves for a short time between two stops.
The boat is stopped if the speed is zero and the next value of stopped is true. This condition captures the cases where the boat stops before moving for some time.
This means that the boat is stopped if the row number is equal to the total number of rows in the group and the speed is zero. This condition captures the cases where the boat stops at the end of the data.
This means that the boat is not stopped if the speed is zero and both the previous and next values of speed are not zero. This condition captures the cases where the boat has a temporary drop in speed but does not stop.
This means that the boat is stopped if the speed is zero and the previous value of stopped is true. This condition captures the cases where the boat continues to stop after stopping for some time.
This means that the boat is stopped if the speed is greater than zero and both the previous and next values of speed are zero. This condition captures the cases where the boat moves for a short time within a stop.
This means that the boat is stopped if the speed is zero and either the previous or next value of stopped is true. This condition corrects any errors like single observations among speed 0 obs groups with stopped=FALSE.
This means that the boat is stopped if the previous value of stopped was false and both the previous and next values of stopped are true. This condition captures the cases where the boat has a false negative in the stopped variable. Another rectifiying measure.
Any NA data in stopped is set as stopped=FALSE.

fishing %>%
  group_by(id,date) %>%
  mutate(seq = ifelse(stopped == TRUE, cumsum(stopped & lag(stopped, default = FALSE) == FALSE), NA)) %>%
  ungroup()->data

data %>%
  group_by(id, date, seq) %>%
  mutate(time_stopped = ifelse(stopped == TRUE, difftime(max(dmy_hm), min(dmy_hm), units = "mins"), NA)) %>%
  ungroup()->data

The data is grouped by id and date to create a variable called “seq” to identify each group of stopped TRUE observations and assign them a serial number from 1 to (n). This is to have a categorical and numerical indicators of each moment when the boats stopped per day.
The data is grouped by id,date and the newly created seq indicator to calculate the time that each vessel remained stopped per day and specific moment. This new variable is called time_stopped

#Data exploration

options(dplyr.summarise.inform = FALSE)


data%>%
  ggplot(aes(time_stopped))+
  geom_histogram()+
  labs(title ="Histogram of time stopped",x="Time stopped",y="Frequency")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 75624 rows containing non-finite outside the scale range
## (`stat_bin()`).

data%>%
  group_by(id,date,seq)%>%
  summarise(ts=mean(time_stopped))%>%
  drop_na()%>%
  ggplot(aes(ts))+
  geom_histogram()+
  labs(x="Mean stopped time (mins)",y="Frequency",title = "Histogram of  mean time stopped")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

data%>%
  group_by(id,Boat.Name,date,seq)%>%
  summarise(ts=mean(time_stopped))%>%
  ungroup()%>%
  ggplot(aes(ts))+
    geom_histogram(bins=50)+
  labs(x="Stopped time (mins)",y="Frequency")+
    facet_wrap(~Boat.Name,scales = "free")

## Warning: Removed 1430 rows containing non-finite outside the scale range
## (`stat_bin()`).

The first graph shows that the highest stopped time is from Bety. The second graph shows a high frequency of 1 minutes stopped.

potential actions: explore extreme data.

data%>%
  filter(!Boat.Name=="Toky")->data
#Filtering out extreme values
#Exploring different time ranges of stopped time


mapbahia+
  geom_point(data=data%>%filter(stopped==TRUE,time_stopped==1),aes(lon,lat),alpha=0.5)+
  scale_fill_distiller(palette="Spectral")

mapbahia+  
  geom_point(data=data%>%filter(stopped==TRUE,time_stopped==2),aes(lon,lat),alpha=0.5)+
  scale_fill_distiller(palette="Spectral")

mapbahia+  
  geom_point(data=data%>%filter(stopped==TRUE,time_stopped==3),aes(lon,lat),alpha=0.5)+
  scale_fill_distiller(palette="Spectral")

mapbahia+
  geom_point(data=data%>%filter(stopped==TRUE,time_stopped>200),aes(lon,lat),alpha=0.5)+
  scale_fill_distiller(palette="Spectral")

The best course of action is to filter out stopped times above 240 mins and equal or less than 1 min.

#Exploring fishing time#

data%>%
  group_by(year(dmy_hm),month(dmy_hm),day(dmy_hm),hour=hour(dmy_hm),Boat.Name)%>%
  summarise(time=mean(time_stopped))%>%
  ggplot(aes(hour,time))+
  geom_smooth()+
  labs(y="Time stopped (mins)",x="Hour of the day (24 hrs)")

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

## Warning: Removed 4288 rows containing non-finite outside the scale range
## (`stat_smooth()`).

data%>%
  ggplot(aes(hour(dmy_hm),time_stopped))+
  geom_smooth(method = 'gam',formula = y ~ s(x, bs = "cc"))+
  facet_wrap(~Home.Location,ncol = 3)+
  labs(y="Average time stopped (mins)",x="Hour of the day (24 hrs)",title = "FIshing time by hour of the day")

## Warning: Removed 59597 rows containing non-finite outside the scale range
## (`stat_smooth()`).

library(leaflet)

## Warning: package 'leaflet' was built under R version 4.3.3

library(leaflet.extras)

## Warning: package 'leaflet.extras' was built under R version 4.3.3

data%>%
  dplyr::select("lng"=lon,lat,time_stopped)->data1

drop_na(data1)->data2

data2%>%
  filter(!lng> -105,time_stopped>10,time_stopped<181)->data2


data2%>%
  group_by(lng,lat)%>%
  summarise(time_stopped=mean(time_stopped))->data2




pal <- colorNumeric("Spectral", domain = data2$time_stopped, na.color = "transparent",reverse = TRUE)



m<-leaflet(options = leafletOptions(maxZoom = 14,minZoom = 9))%>%
  leaflet::setMaxBounds(lng1 = -105.9,lng2 = -105.2,lat1 = 19.8,lat2=21)%>%
  addProviderTiles(providers$Esri.WorldImagery)%>%
  addHeatmap(data = data2,intensity = ~time_stopped,radius = 2,blur=5,max=1,cellSize = 2)%>%
  addLegend(pal= pal, values = data2$time_stopped,
            title="Fishing time (mins)")

## Assuming "lng" and "lat" are longitude and latitude, respectively

Mapfishing

Aldo Zavala

2024-03-27