Mapping Fishing Activities and Suitable Fishing Grounds Using Nighttime Satellite Images and Maximum Entropy Modelling

: Fisheries surveys over broad spatial areas are crucial in deﬁning and delineating appropriate ﬁsheries management areas. Yet accurate mapping and tracking of ﬁshing activities remain largely restricted to developed countries with sufﬁcient resources to use automated identiﬁcation systems and vessel monitoring systems. For many countries, the spatial extent and boundaries of ﬁshing grounds are not completely known. We used satellite images at night to detect ﬁshing grounds in the Philippines for ﬁshing gears that use powerful lights to attract coastal pelagic ﬁshes. We used nightly boat detection data, extracted by U.S. NOAA from the Visible Infrared Imaging Radiometer Suite (VIIRS), for the Philippines from 2012 to 2016, covering 1713 nights, to examine spatio-temporal patterns of ﬁshing activities in the country. Using density-based clustering, we identiﬁed 134 core ﬁshing areas (CFAs) ranging in size from 6 to 23,215 km 2 within the Philippines’ contiguous maritime zone. The CFAs had different seasonal patterns and range of intensities in total light output, possibly reﬂecting differences in multi-gear and multi-species signatures of ﬁshing activities in each ﬁshing ground. Using maximum entropy modeling, we identiﬁed bathymetry and chlorophyll as the main environmental predictors of spatial occurrence of these CFAs when analyzed together, highlighting the multi-gear nature of the CFAs. Applications of the model to speciﬁc CFAs identiﬁed different environmental drivers of ﬁshing distribution, coinciding with known oceanographic associations for a CFA’s dominant target species. This case study highlights nighttime satellite images as a useful source of spatial ﬁshing effort information for ﬁsheries, especially in Southeast Asia.


Introduction
Monitoring and mapping of fishing activities are critical components of planning and management for marine fisheries [1]. Fishing location data have been used to identify and delineate fishing grounds [2], improve stock assessments [1,3], estimate fishing effort [4], and evaluate the impact Anchovies, squids [38] No details available on types of lighting used

VIIRS Boat Detection Data
Nightly VIIRS boat detection (VBD) data from 1 April 2012 to 31 December 2016 were downloaded from the U.S. National Oceanic and Atmospheric Administration National Center for Environmental Information's Earth Observation Group's website (https://www.ngdc.noaa.gov/eog/ viirs/download_phil_boat.html). Individual VBD records represent lit-up pixels in a VIIRS Day/Night Band image [21,30] for one night. The VBD data files contain information on the geo-location of lit-up pixels, radiance value (in nanoWatts/cm 2 /sr), date, and time the image was taken, quality flag, various thresholding values used to distinguish a potential fishing vessel light source versus others, and corresponding VIIRS day/night band image filename. We extracted only VBD points with quality flags 1 (strong boat detection), 2 (weak boat detection), 3 (blurry boat detection), 8 (recurring lights), and 10 (weak and blurry lights), corresponding to radiance spikes that are more likely from marine vessel light sources.
We applied various filtering and aggregation steps to prepare and format the VIIRS boat detection data to the different analyses used to identify core fishing areas, characterize these areas based on cluster analyses, and identify other areas with similar environmental attributes as the core fishing areas in other parts of the country (Figure 1).
Given the satellite's pole-to-pole orbit, some image granules can have zonal overlap especially in mid to higher latitudes with a 2-to 3-hour time difference. To avoid double counting of VBD for the same boat, we identified VBD from two overlapping granules that are 1 km apart and randomly selected one of the two points to represent the vessel's location for that night.
Since most of the fishing gears that use lights (Table 1) catch coastal pelagic fishes [39], we limited our analyses to VBD within the Philippines' 24-nautical mile contiguous zone [40].

Identification of Core Fishing Areas
While fishing activities can occur anywhere in the ocean, previous studies have shown that fishers tend to cluster at specific locations [41,42]. Initial inspection of VIIRS boat detection data revealed similar patterns in the Philippines. We defined these areas as core fishing areas (CFAs) or locations which are repeatedly visited and where there is dense fishing activity.
We accounted for the visitation frequency by creating a standardized 700-m resolution grid and counting the number of years when VBD was detected in each grid cell. We made the resolution slightly smaller than the pixel footprint of the VIIRS satellite images at the nadir to ensure that no grid cell will have more than one VIIRS boat detection for a given night especially near the edges of each image granule. We only used the VBD data within cells that had more than one year of VBD occurrence. This also removed a lot of widely scattered lights. We then obtained the centroid coordinates of each cell with more than one year of VBD data and ran a series of hierarchical density-based spatial clustering of applications with noise (HDBSCAN) [43] using different parameter combinations of 'minimum cluster size' (i.e., 3 to 30) and 'minimum samples' (i.e., 3 to 35), resulting in 1200 different clustering combinations (Supplementary Figure S1). We selected the clustering parameter combination that yielded the least number of points classified as 'noise' while maintaining reasonable clustering of known small pelagic fish fishing areas (e.g., Northeastern Palawan round scad fishing grounds; Tayabas Bay fisheries; and Basilan strait sardine fisheries). Concave hull polygons were estimated for each cluster of points using a Delaunay triangulation algorithm to visualize the clusters. Finally, using high resolution images from Google Earth and shipping traffic data from https://www.marinetraffic.com/ as reference (Supplementary Figure S2), we developed a raster mask to remove polygons in areas with possible high occurrence of non-fishing related shipping activities (e.g., around busy harbors, shipping lanes, and areas with mining activities near the coast) since VBD points in these high-traffic areas could include a large number of non-fishing vessels that turn on additional deck and navigational lights. This mask was also used for the MaxEnt modeling of suitable fishing grounds based on annual VBD occurrence. Annual raster overlays were done using R (v. 3.5.1, R Development Core Team) [44], while hierarchical density-based spatial clustering of applications with noise [45] and concave hull estimations were ran using Python 2.7 (www.python.org).

Clustering CFAs Based on Monthly Patterns and Radiance Values of VBDs
Each CFA can be characterized by the composition of fishing gears operating in the area and the seasonal patterns in fishing. Multiple light-assisted fishing gears (Table 1) most likely operate in each fishing ground. Radiance values can be used as a relative measure of light source intensity with higher radiance corresponding to brighter lights [46]. A narrow radiance range could imply dominance of a single type of fishing gear while a wide range could mean that multiple gears are operating in the area. Seasonal patterns in fishing activities, on the other hand, could be defined by the mean nightly VBD by month, averaged across years within each CFA.
For radiance values, we applied histogram-valued data analysis based on the L2 Wasserstein metric between distributions on log-transformed radiance per CFA [47]. We limited this analysis to VBD points within CFAs occurring during nights with less than 50% moon illumination and less than 50% cloud cover since both factors are known to affect recorded radiance beyond this percentage [21,46] (Supplementary Figure S3). Clustering was performed using the R Package "HistDAWass" [48]. For seasonal patterns in fishing activities within each CFA, we used the mean nightly VBD by month, averaged across years, and standardized the values between 0 and 1, corresponding to the minimum and maximum mean nightly VBD count by month. CFAs were grouped based on hierarchical clustering using Ward's linkage method [49].

MaxEnt Modeling to Identify Environmental Covariates with CFAs
We identified key environmental attributes that define core fishing areas and used these relationships to map out other areas in the Philippines that have similar environmental attributes. Our main assumption is that environmental conditions can be used to partially predict the suitability of a given area for catching small pelagic fishes using light-assisted fishing methods. The relative abundance of target species plays a key role in the selection of fishing grounds by fishers. Spatial distribution of small pelagic fishes has been known to correlate with various environmental attributes that link with food availability and habitat suitability. We recognize, however, that other socio-economic (e.g., cost of fishing, distance from ports), cultural (e.g., adherence to traditions and/or familial knowledge), governance (e.g., access and restrictions to certain areas), and safety factors are also considered by fishers in choosing their fishing grounds [50].
We applied MaxEnt version 3.4.1 [51], a machine learning algorithm based on the concept of maximum entropy [52], to determine the environmental conditions that characterize core fishing areas and identify other areas in the country with similar suitable environmental conditions that could also have high fish productivity. Since the VIIRS satellite only takes one image per night and small pelagic fishes are also caught using other gears that do not employ powerful lights, we conservatively treated the VBD locations as presence-only data and randomly derived pseudo-absences or background locations outside the CFAs. MaxEnt is one of the few models available that can handle presence-only data and has been used extensively for species distribution modeling [53,54]. It estimates the relative suitability of one place versus another by comparing the conditional density of covariates at presence sites with the unconditional density of covariates across the study area [55]. The probability of occurrences are returned which is relative to the 'typical' conditions generated from available presence points.
We applied MaxEnt at an aggregated 4 km resolution to match the resolution of remotely-sensed environmental datasets used as predictors ( Table 2). The model domain is the Philippines' 24 nautical mile contiguous zone bounded from 4 • to 22 • North and from 116 • to 128 • East. 'Presence' locations were defined as the centroid coordinates of 4-km grid cells that had VBD data and were inside any of the CFAs for the periods 2013 to 2016. VBD data for 2012 were incomplete (i.e., starting only in April 2012) so they were not used in the MaxEnt models. Background points (~10,000) were randomly selected within the model domain. MaxEnt's regularization parameter, controlling for over-fitting, and the selection of feature classes (e.g., linear, quadratic, product, threshold, or hinge) [55] were tuned for each model using the ENMevaluate function in the ENMeval package in R [56]. For each environmental predictor, the 'auto feature' selection method was used to allow MaxEnt to determine the appropriate feature class to use for each environmental predictor in each model with 'maximum iterations' set at 5000.
We collated environmental covariates commonly associated with the distribution of small pelagic fishes [57][58][59][60] (Table 2; Supplementary Figure S4). Gridded bathymetry was downloaded from the General Bathymetric Chart of the Oceans (GEBCO) website. Monthly netCDF files for the other environmental predictors were downloaded from the NOAA Coast Watch's Environmental Research Division's Data Access Program archives for the years 2013 to 2016. Climatological and annual means were calculated for each predictor across the four years. All raster layers were projected to a Behrmann's equal area projection at 4km resolution using bilinear interpolation under the 'raster' package in R [61].

MaxEnt Scenarios
Using the biomod2 package in R [62], we ran three different sets of MaxEnt models and scenarios: (a) climatology, (b) annual, and (c) specific CFAs (Table 3). For the climatology runs, we generated time-series averaged rasters for all six dynamic predictors (i.e., except bathymetry) from 2013 to 2016. Under this set, models were run with and without bathymetry as a predictor and accounting for possible effects of monsoons which have been known to affect small pelagic fish production in the Philippines [63].
The exact location of fishing activities, even within the same CFA, often changes across years. To account for annual variability of dynamic predictors on annual aggregated locations of VBD inside CFAs, we developed annual MaxEnt models from 2013 to 2016, using respective annual VBD data and annual averages for each predictor. Since the CFAs most likely differ in terms of the fishing fleet composition and targeted fish species, we selected two CFAs for which the dominant target species group is known and where most fishing gears that use lights are almost exclusively targeting these species. These CFAs are (a) the round scad fishery in Northeastern Palawan island (CFA # 106 in Figure  3) and (b) the sardine fisheries in the Sulu Archipelago (CFA # 70).
To identify other areas within the Philippines' 24 nautical mile contiguous zone with environmental characteristics similar to the delineated CFAs, we generated predicted CFA suitability maps for all 10 MaxEnt models (Table 3). Mean predicted suitability values per model was generated from suitability maps produced in each of the cross-validation runs (n = 10 per model). The resulting suitability distributions were compared with the location of existing CFAs.  Figure 3) and climatology of environmental predictors C2. Sulu VBD in CFA # 70 (Sulu; Figure 3) and climatology of environmental predictors

Model Validation and Evaluation
To test each MaxEnt model's predictive accuracy, we applied two repetitions of a 5-fold cross-validation that splits the presence data into five subsets. For each run, four of the five subsets (i.e., 80% of presence data) were aggregated and used to train the model, the remaining subset (i.e., 20% of presence data) was used to evaluate how well the trained model can accurately predict these occurrences [64]. Results were averaged across the repeated 5-fold cross-validation runs to predict the current distribution of suitable fishing grounds for light-assisted fisheries within the Philippines' 24-nautical mile contiguous zone. Accuracy was measured for each cross-validation run using the area under the curve of the receiver operating characteristic (AUC) and the true skill statistic (TSS) [65][66][67]. AUC and TSS values range from 0 to 1. AUC values greater than 0.5 indicate that model predictions are better than random predictions. TSS values greater than 0.4 indicate good models [68].
We reported 'variable importance' as a measure of the contribution of variables to the discrimination of presence and pseudo-absence points. We also ran jackknife tests of environmental importance (i.e., leave-one-out) which identifies improvements in model fit by adding or removing each individual variable and measuring changes in training or testing data gain. Response curves were generated and examined to derive the preferred range of environmental variables of potential fishing areas.

Results
VIIRS boat detection (VBD) locations within the Philippine's 24-nautical mile contiguous zone are highly clustered with specific areas being frequently visited (Figure 2a). Spatial patterns are also apparent in terms of mean radiance values with most VBD locations having radiance values greater than 2.0 nanoW/cm 2 /sr (Figure 2b). VBD with less than 2.0 radiance values are clearly clustered together in some areas. Across years, the spatial patterns of VBD aggregations within the Philippines has not change substantially from 2012 to 2016.

Core Fishing Areas in the Philippines Maritime Contiguous Zone
A total of 694,793 VBD points were recorded within the Philippines' contiguous zone over 1713 nights. Eighty-five percent (85%) of these data points were in cells that were visited for at least two years over the five-year period (Figure 3).
We identified 134 Core Fishing Areas (CFAs) for light-assisted fishing in the Philippines based on HDBSCAN results using a 'minimum cluster size' of 4 and 'minimum samples' value of 25 (Figure 3; Supplementary Table S1). This HDBSCAN parameterization initially labelled 94% of the frequently visited VBD cell centroids into 160 clusters. We removed clusters that were (a) too small, (b) included or were located near major ports / piers, and/or (c) located in major shipping lanes.
The 134 CFAs ranged in size from 6 km 2 to 23,215 km 2 with a mean size of 700 km 2 . The maximum number of VBD for a given night was detected in March 3, 2014 with 1,124 lit-up pixels within the entire Philippines Contiguous Zone, 87% or 974 of these pixels were within CFAs. Per CFA, the maximum number of lights detected for a single night ranged from 2 to 227. Figure 3. Core fishing areas (n = 134) identified from HDBSCAN and after removing clusters in areas that overlap with significant vessel traffic and close to major ports and piers, including mining jetties. Details for each CFAs and corresponding names can be found in Supplementary Table S1. Seasonal patterns in fishing activity data varied per CFA (Figure 4 and Supplementary Figure S5-a). Cluster analysis applied to averaged monthly VBDs showed two different groups based on the months when peak fishing activity occur, generally corresponding with the monsoon months: (a) northeast monsoon (November to April) and (b) southwest monsoon (May to October).

VBD and Environmental Variables
Mean sea surface temperature, chlorophyll a, and bathymetry inside the CFAs differed from those outside ( Figure 6). CFAs have mean depth ranging between 5 m and 1,480 m, mean sea surface temperature between 28.5 • C and 30.6 • C, and mean surface chlorophyll a values between 0.2 mg/m 3 and 2.2 mg/m 3 . The difference between monsoons was apparent only in mean chlorophyll concentrations with the VBD points during the southwest monsoon located in areas with higher mean chlorophyll a concentrations (0.75 mg/m 3 for southwest monsoon vs. 0.43 mg/m 3 for northeast monsoon; Welch two sample t-test p < 0.01). All MaxEnt models performed better than random (i.e., AUC > 0.5 and TSS > 0.4; Table 4). Among all the models, the one without bathymetry had the lowest AUC and TSS scores. The full climatology, monsoonal, and annual models all had closely similar accuracy scores. The CFA-specific models performed much better than the other models, with high AUC and TSS scores (greater than 0.9 for both), highlighting the unique environmental attributes of these CFAs. Bathymetry ranked highest in variable importance among all non-CFA-specific models, with permutation importance values ranging from 62% (2015) to 87% (Full model) (Figure 7). This was followed by chlorophyll a with permutation importance ranging from 2% (Southwest monsoon) to 20% (2015). In the model without bathymetry, chlorophyll explained most of the difference between CFA and non-CFA locations with a mean permutation importance of 54%. Jackknife tests also identified bathymetry as the variable with the highest model training, testing, and AUC gain when used in isolation (i.e., most useful information by itself) and the greatest decline in gain when removed from the model (i.e., contains most information that is not present in the other variables). Sea surface salinity and significant wave height ranked third in terms of variable permutation importance, particularly in the monsoonal models. Sea surface temperature and sea surface height varied in their contributions to the model but mostly at less than 10% importance. In contrast to this general pattern, the variables contributing to the CFA-specific models were more even with bathymetry ranking third from the bottom in importance for the sardine fishery in Sulu (CFA # 70). In the Northeast Palawan model, sea surface salinity ranked highest in importance (42%) followed by bathymetry (29%) and chlorophyll (11%). In the Sulu CFA model, chlorophyll a ranked highest in permutation importance (73%) with the other variables each having less than 10% permutation importance. Based on the response curves from the MaxEnt full climatology model, the logistic output or probability of occurrence of a site becoming a core fishing area increased with depth from the shallowest depths up to around 100 m then decreased rapidly to depths of 500 m and gradually declined further with depths deeper than 4000 m. The probability also peaked with chlorophyll a values between 0.3 and 2.0 mg/cm 3 . Frequently visited fishing sites also occurred more in areas with sea surface temperatures between 28.6 • C and 30.6 • C. The effect of sea surface salinity, sea surface height, and mixed layer thickness in discriminating CFA versus non-CFA areas was smaller compared to the first three variables mentioned as indicated by the almost flat response across their entire range of values.
The predicted suitability of marine areas for light-assisted fishing based on the climatology and annual MaxEnt models are shown in Figure 8. Predicted suitable areas with environmental conditions typical of the CFAs were found in shallow, nearshore areas around existing CFAs. Predicted suitable areas did not vary greatly across monsoons and years. For the CFA-specific MaxEnt models, we found highly constrained areas of high suitability. However, the predicted areas corresponded with other CFAs that were not included in the model (e.g., CFA # 6, 18, 80-94, and 102 for the Northeast Palawan round scad CFA and CFA # 39, 53, 59, 67, 69, 71, and 72 for the Sulu sardine CFA).

Discussion
The VIIRS nighttime imaging capabilities significantly enhanced those of the DMSP satellites, particularly in terms of radiometric and spatial resolution, allowing for the nightly extraction of locations and radiance information of maritime lights [21,32]. We applied this freely accessible and novel dataset to identify, delineate, and characterize core fishing areas (CFAs) and determine the suitability of other parts of Philippine coastal waters to light-assisted fisheries and their target species.
Many of the large CFAs identified are known primary fishing grounds for various coastal pelagic fishes. The Sulu CFA (#70) in Southern Mindanao is dominated by sardine purse seine fishing. It is one of the largest sources of the country's sardine (Sardinella lemuru) supply and has been subjected to seasonal fishery closures between November and March since 2011 due to concerns about overfishing and stock declines [69]. The Northeast Palawan CFA (#106) and those west of Palawan (e.g., CFAs # 102 to 117) are known round scad (Decapterus macrosoma and Decapterus russelli) fishing areas, also subjected to annual seasonal closure since 2015 [70]. Visayan Sea 1 (#37) and 2 (#40) CFAs are also part of a small pelagic fish seasonal closure established back in the 1970s primarily for sardines (Sardinella gibbosa).
A similar study was done globally, identifying clusters or 'agglomeration areas' of marine lights from cloud-free monthly composites of the VIIRS day/night band images [1]. The authors limited the presentation of their results to only the top 100 agglomeration areas by size. Within this list, they identified six (6) agglomeration areas for the Philippines. All these six agglomeration areas were also identified as CFAs in our study. Twenty-one of the 134 CFAs are larger than 1000 km 2 . Unfortunately, information on the Philippine agglomeration areas were not available in Zhao et al. [71] since they only presented statistics for the top 50 of their top 100 agglomeration areas. The difference in the number of agglomeration areas and CFAs identified for the Philippines could be due to Zhao et al.'s [71] use of only cloud-free monthly composites and a fixed distance threshold (i.e., 10 km) to group pixels.
Bathymetry, chlorophyll a, and sea surface temperature have all been cited as predictors of small pelagic fish spatial distribution in previous studies. Catches and fishing grounds for the Indian mackerel (Rastrelliger kanagurta) in Indonesia correlated significantly with chlorophyll a and sea surface temperature [72]. Bottom depth and sea surface temperature explained variability in European sardine (Sardinella pilchardus) presence in the Mediterranean Sea while bottom depth and surface chlorophyll concentration predicted European anchovy (Engraulis encrasicolus) distribution [73][74][75]. These three variables, along with salinity, are also used in mapping global distributions of fish species through Aquamaps (www.aquamaps.org) [76,77]. These variables influence suitable habitats for small pelagic fishes through zooplankton production and availability [78].
Bathymetry consistently ranked highest in discriminating between CFA and non-CFA sites, accounting for the most variability not captured in the other variables. Chlorophyll accounted for most of the variability in the MaxEnt model that excluded bathymetry indicating some correlation between chlorophyll a and bathymetry (Pearson's correlation r = 0.25; p < 0.01). Indeed, chlorophyll a tended to decline with increasing depth. However, zooming in on some of the CFAs shows how VBD data track, with surprising precision, bathymetry, and depth contours using a 1-km resolution bathymetry data (Figure 9). Removing bathymetry as a predictor resulted in a lower accuracy model (AUC from 0.88 to 0.84 and TSS from 0.59 to 0.52; Table 4) and broader predicted areas of suitability ( Figure 8). Aggregating data at a 4-km resolution removes these fine-scale spatial patterns in nighttime lights and fishing sites. Finer resolution environmental data could explain more of the variability in the distribution of nighttime lights in CFAs. Unfortunately, current available monthly remote sensing data for chlorophyll and sea surface temperatures for this region are limited to a 4-km resolution (i.e., MODIS images) while finer oceanographic models are not available for the entire study area. The dominance of bathymetry in distinguishing CFA from non-CFA sites could be attributed to (1) the depth preference of target species and (2) the multi-species nature of fisheries within CFAs. A search in Fishbase (www.fishbase.org) on depth information for the dominant small pelagic fishes caught in Philippines waters shows that commonly targeted small pelagic fishes by the Philippine fishing fleet are often found in depths of less than 200 m [63]. In addition, CFAs most likely have a large mix of different fishing gears that target various species of small pelagic fishes. The presence of lights in deeper areas, although not dense enough or frequently visited to classify as a core fishing area, indicate that fishing vessels using lights are able to fish further offshore and not constrained by vessel capability. Since different species or species group tend to respond differently to water column characteristics (e.g., SST and chlorophyll i concentrations), the common denominator in their environmental preference becomes bottom depth, a static variable which then gets captured as the main predictor of CFA versus non-CFA sites. If the CFAs could be differentiated according to their dominant target species, the effect of dynamic predictors in differentiating CFAs from non-CFAs for that particular species could be bigger than bathymetry.
One of the main uses of the identified core fishing areas and habitat suitability map in this study is to delineate fisheries management areas as part of implementing ecosystem approaches to fisheries management (EAFM) [8,9]. While we have provided here a way to identify frequently visited fishing grounds and predict habitat suitability of non-fishing grounds based on nighttime lights activity and environmental attributes, delineation of meaningful fisheries management unit as part EAFM should be done in consultation with relevant stakeholders [2]. Some of the CFAs could be grouped into larger fisheries management areas based on similarities in their target species and fishing gears used. In the absence of actual fishing data available for each CFA, the seasonality of fishing activities and the distribution of VBD radiance values could be used as a proxy for these characteristics. For example, the northern Palawan CFAs (e.g., #102 to 117) can be grouped into one FMA since they share the same seasonality in fishing activities (during northeast monsoon; Figure 4) and radiance signature ( Figure 5; see also Supplementary Figure S5 for maps of the CFAs color-coded according to cluster membership). The CFAs west of Mindoro island (#s 122 to 125) could also be considered as a separate fisheries management area given their similarities in both seasonality of fishing activities (during northeast monsoon) and radiance signature (mostly low radiance values). Of course, other fishing activities and practices should also be considered in defining appropriate fisheries management areas since the VBD only covers fishing vessels that use lights to attract fish. As Jennings and Lee [2] stated, the ultimate delineation and definition of 'fishing grounds' would be a societal decision of choosing cut-offs and satisfying potential trade-offs in the use of limited marine space.
The VBD data could also provide important information on the levels of fishing effort over time. The estimated number of fishing vessels using lights in the Philippines from the VBD data (~1,000) is most likely an underestimate since the VIIRS instrument takes only one image per night and for most nights, there is always a part of the country that is obscured by thick cloud cover. Summing up the maximum nightly VBD per CFA for selected larger, well-separated CFAs (i.e., #s 19,66,70,75,106,122) already results in 925 estimated light boats or fishing vessels with on-board lights. In comparison, there are an estimated 6,371 commercial fishing vessels (i.e., 3 gross tons and heavier) in the country as of 2007 [79]. Despite this potential underestimation of absolute fishing effort, the quarterly pattern of nighttime fishing lights still tracks closely the national statistics of landed catches of small pelagic fish (Supplementary Figure S6). Monthly VBD data show a secondary peak between the third and fourth quarter of each year which is not captured in the quarterly landed catch data. This secondary peak in the VBD data is most likely associated with the varying seasonal patterns of fishing activities among CFAs with some CFAs' peak VBD counts occurring around this period ( Figure 4). Subject to further ground validation, the temporal trends in VBD data could be used as a proxy for relative fishing effort at the national level for small pelagic fish catch, a key input for normalizing catch data.

Limitations and Future Improvements
In using the VBD data for fisheries studies, additional filtering is needed to ensure that most of the data are fishing vessels. Although VBD use quality flags to identify non-vessel light sources in the data, it has no flags to distinguish fishing from non-fishing vessels (e.g., cargo ships, tankers, passenger ships, construction, and mining-related temporary structures, among others). Navigation routes are clearly seen from VBD overlays which correspond with known shipping lanes (e.g., between Mindoro and Panay Islands; Sabah, Malaysia and Tawi-Tawi, Philippines; Figure 3 and Supplementary Figure S1). Major ports and mining areas also have high density aggregations of nighttime lights. It is possible that many of these lights are also fishing vessels but without ground validation, we decided to exclude them as core fishing areas and in the MaxEnt modeling of suitable areas. Even though we removed some of these potential CFAs, the fishing ground suitability map produced from MaxEnt was still able to identify some of them as potential fishing grounds for small pelagic fishes based on similarities in typical ocean conditions to most of the CFAs where lights are found (Figure 3; e.g., northeast tip of Mindanao).
We ran the maximum entropy model under the assumption that areas where fishers often fish (e.g., CFAs) are also areas of high relative abundance for the target species group, in this case, small pelagic fishes, and that environmental factors can be used to predict habitat suitability [80]. At a resolution of 742 m, it is surprising that many of the cells where nighttime lights are found in the last five years are revisited over the years. We hypothesized that environmental factors played a large role in the spatial structure of the nighttime lights, particularly for these frequently-visited and high-density sites. However, other factors also influence fishers' decisions in selecting fishing grounds such as projected fishing costs, estimated risk, and adherence to traditional fishing grounds [50]. These are not captured in our model.
The utility of the VBD data in informing ecosystem approaches to fisheries management could be enhanced by studies that can identify the relationships between radiance values and fishing characteristics (e.g., fishing gear or strategy). This would help translate the VBD products from remote sensing outputs to the more common fisheries metrics (e.g., number of vessels, fishing gears, and/or species). In addition, the effect of cloud cover on detection probability and limits should be investigated further especially when utilizing the VBD data for near real-time applications. By aggregating the boat detection data annually or over several years, we minimized the potential spatial bias brought about by spatial differences in the frequency of cloudy nights. Supplementary Figure S7 in the supplementary file shows that lights are still detected even in areas with high frequencies of cloudy nights. Despite these limitations, having nightly information on distributions of light-assisted fishing activities freely accessible and available to fisheries managers partially fills the crucial spatial data gap for implementing ecosystem approaches to fisheries management.

Conclusions
Any marine resource management requires an understanding of the spatial and temporal patterns of both resource use and target resources and the underlying processes. Yet collecting information on spatial and temporal patterns of resource use at relevant scales for locally-adapted action remains sparse, especially in marine fisheries in developing countries. We used the VIIRS boat detection data service from NOAA to identify high density fishing grounds and predict areas within the Philippines's contiguous zone with similar environmental conditions in aid of implementing ecosystem approaches to fisheries management. Nighttime lights are extensively detected in the Philippine coastal waters. These lights spatially cluster and we have identified 134 distinct core fishing areas. Aggregating all boat detections inside CFAs, these areas are strongly distinguished from non-CFA areas mainly by bathymetry, owing most probably to the multi-gear nature of fishing and the varying fishing operation strategies employed within each CFA. In two of the largest CFAs which we know are each targeting round scads and sardines, variable importance in the MaxEnt models corresponded with environmental parameters previously known to be associated with these species' occurrence. This highlights the potential model improvement if the catch composition and fishing fleets within each CFA can be characterized. In addition, using finer spatial resolution environmental predictors could enhance model precision and differentiate suitable fishing areas even within the same CFA.
The VIIRS boat detection data is provided freely and for public use. Its applications can be further expanded by conducting ground-truthing and calibration studies to relate VBD radiance values to common fishery statistics (e.g., number of vessels, types of vessels, comparative fishing effort based on boat size, etc.). However, even with the current dataset, having nightly coverage of fishing activities, even if just a subset of the entire fishing sector, provides an invaluable data stream of fishing effort for data-deficient fisheries, especially in Southeast Asia. Our results could also be used as a basis for establishing fisheries management areas in the Philippines and our approach could be applied to other countries that have high numbers of light-employing fishing vessels.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Figure S1: HDBSCAN parameter combination; Figure S2: Marine traffic density map; Figure S3: Mean radiance versus percentage cloud cover and moon illumination; Figure S4: Mean climatology maps of environmental predictors; Figure S5: Core fishing areas mapped according to seasonal and radiance clusters; Figure S6: Comparison of VBD trend and reported quarterly production of small pelagic fish in the Philippines; Figure S7. Spatial frequency of cloud cover in relation to VIIRS boat detection data; Table S1: Characteristics of core fishing areas.