Update (2019-09-21): I revisited this post and I tried to redo the plot but instead I used base R graphics plus the mapplots package and the data.table package for manipulation.

Last May 24-25, 2018, I attended the 9th West Pacific East Asia - National Stock Assessment Program (WPEA-NSAP) Tuna Catch Estimates Review Workshop at Ocean 101, Catangnan, General Luna (Siargao), Surigao Del Norte. The said activity aimed to review the NSAP port sampling data collected in each region and what information the BFAR regional offices had compiled that could be used in the annual catch estimates of tuna species. One of the presentations from another region interests me - they showed a map of the fishing grounds in their region and overlaid a pie chart showing the landed catch composition of tuna. When I returned to the office, I tried to make the same graph but using the R language.

Note: Since we cannot share the NSAP data without the permission of the higher officials, I used the data from the Philippine Statistics Authority (PSA) and is freely accessible here. I pre-formatted the .csv file and you can download it here.

The data

As noted above, the data I used was extracted from PSA and I only selected the 2015 to 2017 data of three tuna species - bigeye tuna, skipjack, and yellowfin tuna. Furthermore, I only processed the tuna production of marine municipal fisheries. Lastly, only the data from the six provinces of Bicol Region was used - Albay, Camarines Norte, Camarines Sur, Catanduanes, Masbate, and Sorsogon.

The process

To begin, let us load the needed packages.

library(mapdata)
library(mapplots)
library(data.table)

Data manipulation

Before going to making the graph, let us first look at our data.

municipal <- fread("data/raw-data/municipal.csv", header = TRUE)

head(municipal)

Now we will add a new column containing the sum of columns 2015, 2016, and 2017. We will be going to sum it in a row-wise fashion.

municipal[, total := rowSums(.SD), .SDcols = 3:5]

head(municipal)
##           province        species    2015    2016    2017   total
## 1:           Albay    Bigeye tuna  285.74  281.94  334.97  902.65
## 2:           Albay       Skipjack    3.24   54.83  114.55  172.62
## 3:           Albay Yellowfin tuna  262.57  262.66  289.43  814.66
## 4: Camarines Norte    Bigeye tuna  448.39  627.51  620.60 1696.50
## 5: Camarines Norte       Skipjack 1347.72 1381.78 1215.91 3945.41
## 6: Camarines Norte Yellowfin tuna  713.83  474.89  382.25 1570.97

Now, we will only select the province, species, and the total columns.

municipal <- municipal[, .(province, species, total)]

We will convert the data into a wide format to add the coordinates and convert again to long format (I still don’t know the easiest method to do this).

# Convert to wide format
municipal_wide <- dcast(municipal, province ~ species, value.var = "total")

# Add the coordinates of each province
municipal_wide[, `:=` (lon = c(123.632377, 122.763304, 123.348615, 124.242160, 123.558856, 123.930399), 
                       lat = c(13.171977, 14.139026, 13.525020, 13.708868, 12.306024, 12.759986))]

# Then convert back to long format
municipal <- melt(municipal_wide, id.vars = c("province", "lon", "lat"))

head(municipal)
##           province      lon      lat    variable   value
## 1:           Albay 123.6324 13.17198 Bigeye tuna  902.65
## 2: Camarines Norte 122.7633 14.13903 Bigeye tuna 1696.50
## 3:   Camarines Sur 123.3486 13.52502 Bigeye tuna  397.27
## 4:     Catanduanes 124.2422 13.70887 Bigeye tuna  716.94
## 5:         Masbate 123.5589 12.30602 Bigeye tuna 1314.28
## 6:        Sorsogon 123.9304 12.75999 Bigeye tuna  362.89

The pie charts above the map

We will use the function draw.pie from the mapplots package. The documentation for this can be found at https://rdrr.io/cran/mapplots/man/draw.pie.html.

The final output is shown below:

Using the mapplots::draw.pie function to add pie charts on map.

As we can see, more tuna are landed in Camarines Norte compared to other provinces and the dominant species is the skipjack (Katsuwonus pelamis). In Camarines Sur and Catanduanes, the dominant species is the yellowfin tuna (Thunnus albacares), whereas bigeye tuna (Thunnus obesus) dominates the landed catch in Masbate.

That’s it.

The code to generate the plot

# The area of the Bicol Region;
xlim <- c(122, 125)
ylim <- c(11.5, 14.5)

# Creates an xyz object for use with the function draw.pie
xyz <- make.xyz(municipal$lon, municipal$lat, municipal$value, municipal$variable)

# Colors used
col <- c("#003366", "#CCCC66", "#CC3366")

# The plot of the pie chart above the map
tiff("pie-on-map.tiff", width = 8, height = 5.5, units = "in",
     res = 200, type = "cairo")
par(mai = c(0.5, 0.5, 0.35, 0.2), omi = c(0.25, 0.5, 0, 0),
    mgp = c(2.5, 0.5, 0), family = "Liberation Serif")
basemap(xlim = c(121, 126), ylim = c(11.5, 14.5), bg = "white",
        main = "Distribution of three major species of tuna in Bicol Region")
map('world2Hires', xlim = xlim, ylim = ylim, add = TRUE)
draw.pie(xyz$x, xyz$y, xyz$z, radius = 0.3, col = col)
legend.pie(121.5, 11.75, labels = c("Bigeye", "Skipjack", "Yellowfin"), 
           radius = 0.2, bty = "n", col = col, cex = 0.8, label.dist = 1.3)
text(121.5, 12.1, "Tuna Species:", cex = 0.8, font = 2)
text(123.4, 14.3, "Camarines\nNorte", cex = 0.8, font = 2)
text(124.8, 13.75, "Catanduanes", cex = 0.8, font = 2)
text(122.9, 13.5, "Camarines\nSur", cex = 0.8, font = 2)
text(124.1, 13.25, "Albay", cex = 0.8, font = 2)
text(124.3, 12.8, "Sorsogon", cex = 0.8, font = 2)
text(123.5, 12, "Masbate", cex = 0.8, font = 2)
mtext("Data Source: Philippine Statistics Authority", side = 1, outer = TRUE,
      adj = 1, cex = 0.8, font = 3)
dev.off()