Project functions: Cloud mask for Landsat 8: maskL8sr() Cloud mask for Sentinel 2: maskS2clouds() Add EOM1 index to Sentinel 2: addEOMI1() Add EWI index to Sentinel 2: addEWI() Add NDWI index to Sentinel 2: addNDWI() Add EWI index to Landsat 8: l8_addEWI() Add NDWI index to Landsat 8: l8_addNDWI() Add NWI index to Landsat 8: l8_addNWI() Add 5:4 band ratio to Landsat 8: l8_5_4ratio()

Yucatán and CAFO lagoons

Sentinel-2 2019 composite image of Yucatán. Turquoise polygons represent CAFO Lagoons. S2 image from Google Earth Engine. Lagoons from story map.

1. Introduction (CB)

Concentrated Animal Feeding Operations, or CAFOS, are intensive farming operations in which large numbers of pigs (in some cases up to 80,000 pigs) are held in indoor pens for meat production. Large amounts of fecal waste are generated from these animals and are disposed of in “lagunas,” or open pools of water and pig feces. These pens are breeding grounds for dangerous bacteria like Salmonella and E. coli and are not only unsanitary environments for the animals, but for communities living nearby.

Industrial pig farming has been growing in Mexico’s Yucatán Peninsula with adverse environmental and social impacts, including; the destruction of livelihoods (mainly due to environmental degradation), increased rates of asthma, cancer, premature births, low birth weights, deforestation, and overconsumption and contamination of water. As part of our final project for Geospatial Analysis with R, we will be building upon previous work by Karen Vazques Hudlet, Michael Cecil, and Manuel Llano and focused on studying the growth of Concentrated Animal Feeding Operations (CAFOs) in the Yucatán peninsula using Rgee for RStudio and Landsat 8 and Sentinel 2 images from Google Earth Engine. We developed a method for identifying CAFO’s refuse lagoons from satellite imagery— as these lagoons have a unique spectral signature compared to other bodies of water— to identify when the lagoons first began to appear in the peninsula and at what scale and rate their number has grown over time. We used both Sentinels 2 and Landsat 8 satellite images in our analyses. We created an image collection for each year that included every image from the span of a few months (one image for each day) and used a cloud mask to remove cloud interference. We trained a Random Forest classifier for both Sentinel 2 images collections (2019-2022) and Landsat8(2014-2018) using spectral signatures from existing lagunas and sample polygons of “not lagunas” (polygons of mixed spectral signatures including, forest, urban, beach, etc.). We used our trained classifier to loop through each year to find and classify lagunas and represent the data in R with maps and graphs that show the distribution of CAFOS in the Yucatan Peninsula. We hope this information will aid in identifying and mapping CAFOS in Yucatan even if they are previously undocumented and will add to previous research on CAFOs in Yucatán and support community claims of increasing, illegal, industrial swine operations in their territories.

2. Methods (AD)

The methods used in this project are recorded in the project_work.Rmd file in the vignette folder. We refer to the chunks in that file to describe our work flow.

This project uses the rgee package in Rstudio to access Google Earth Engine. Most of the processing was done on Earth Engine servers. While some coding was done on the Earth Engine code editor, most of it was done in Rstudio.

Our workflow was the following:

Uploaded study area shapefile to use in earth engine. We used a state polygon layer from a Mexican government website. (chunk 1)
Uploaded shapefile of CAFO lagoons to use in earth engine. The polygon layer was provided by Mike Cecil from his previous research for the story map project. (Chunk 2)

Project study area: The state of Yucatán and identified CAFO lagoons.

Used rgee to extract a cloudless satellite image composite of the state of Yucatán. We used Sentinel 2 because of the higher 10 meter resolution. We chose images from April because it is the driest month in the region. The images are from 2019 because that is the latest lagoon date recorded in the lagoon shapefile. We sourced cloud mask functions for removing cloud reflectance values in both Sentinel 2 and Landsat 8 images. The cloud mask functions are saved as part of the cafos23 package.The median cell value was then taken from every image included in the time frame (April 10 - 20) to produce a cloudless image. (chunk 3, see functions in R folder)

Dates recorded for CAFO lagoons included in the CAFOS shapefile. These years are based off of land cover change and are not all accurate, especially for the 2000 and earlier lagoons. This could explain why the data is intensely skewed.

We used a true and false color composite (NIR, Red, Green) of the 2019 Sentinel 2 image to examine which lagoon polygons were actual CAFO lagoons in 2019. The FCC was the most useful for identifying lagoons since they appeared black and in dark contrast to the surrounding vegetation, which appeared red. While most of the polygons were indeed lagoons, we chose to only use lagoons that were not too dry, covered by too much vegetation, or if more than 40% of the polygon was not covering the lagoon. This helped ensure a unique spectral signature for lagoons that would help differentiate it from bare soil and vegetation. We rejected 85 out of 320 lagoons, keeping 235 lagoons that we refer to as known 2019 lagoons. (chunks 4, 6, 7)

CAFO Lagoon (center, index 35) over 2019 S-2 FCC image. This lagoon was included in the 2019 lagoon count.

CAFO Lagoon (center, index 273) over 2019 S-2 FCC image. This lagoon was rejected from the 2019 lagoon count since it looks to be mostly vegetation and soil

After identifying which lagoons to reject, we created a new lagoon feature collection (lagunas2019_fc) that had only the known 2019 lagoons. (chunks 8, 9)
We then took random samples from the known 2019 lagoons and areas that were not lagoons in 2019. The non-lagoon sample points came from manually created polygons on the Earth Engine code editor. The Sentinel-2 2019 composite image was sampled. The sample points were merged together, each point had one attribute (property) labeled as presence. Presence 0 meant the point was not from a lagoon, and presence 1 meant the point was from a lagoon. (chunks 10, 12, 13)

Sample points in red used for training Random Forest model. Turquoise polygons are lagoons, yellow polygons are non lagoons.

The sample points were used to extract the reflectance values from the Sentinel-2 2019 image of the lagoons and non-lagoons. We used only the values from bands 4 (Red) and 8 (NIR). Additionally we sourced functions to add three indices to the points:

NDWI (Normalized Difference Water Index: Derived from NIR and SWIR and sensitive to changes in liquid water content and in spongy mesophyll of vegetation canopies)
EWI (Simple Enhanced Water Index: (Green-(NIR+MIR))/(Green+(NIR+MIR))
EOMI1 (Exogenous organic matter index: Developed for detecting manure application on fields)

The bands’ and indexes’ values were used to train an 1000-branch Random Forest model to create a classifier for the spectral signatures of the lagoon and non-lagoon classes. (chunk 11, 14, see functions in R folder)

The trained classifier was then used to classify the 2019 Sentinel-2 image into either lagoons or non-lagoons.

2019 classification with Sentinel-2 trained classifier. Blue areas are classified as lagoons. Tan areas are classified as non-lagoons.

We then overlayed the original 320 lagoons over the 2019 classification. We used the ee.function reduceRegions with a mean reducer to calculate the average cell value of the pixels within each lagoon polygon. The cells can either have a value of 0 (non-lagoon) or 1 (lagoon), so the closer the average is to 1, the greater the proportion of the polygon that is classified as lagoon. The image below demonstrates a lagoon polygon that is 58% classified as lagoon. (chunk 16)

58% of the red polygon (index 315) is classified as lagoon. When viewed on the color composite, there seems to be a lot of vegetation obscuring the lagoon, making it unclear if it is indeed a lagoon. This lagoon would be rejected by the classifier cut-off proportion that we chose (75%).

We tested various proportions to be used as the cut-off for which lagoon polygons should be classified as lagoon, before settling on 0.75. We found that if we classified only lagoon polygons that had at least 75% of their areas classified as lagoons (blue), this maximized the number of known lagoons that we correctly classified and minimized how much we overclassified. With a 75% cut-off, we classified 248 of the original 320 lagoon polygons as lagoons. However, this was 13 additional lagoons to the 235 known 2019 lagoons. We cross-checked the index of the lagoons with the list of lagoons to be removed from 2019 and found that we misclassified 21 non-lagoons as lagoons. The histogram below displays the distribution of the lagoon polygon areas’ classifications. The S-2 classifier has the tendency to overclassify the lagoons, which is why we chose to go with such a high proportion cut-off of 75%. (chunk 17)

At this point we had planned to classify and calculate the lagoon count of older Sentinel-2 Yucatán images, however we could not find good images of the region in the ee collection, which only went back to 2018 in Yucatán. So instead we decided to use the Landsat 8 collection, which had good images of Yucatán back until 2014.
We repeated the previous steps to train a Random Tree classifier with a Landsat 8 composite image from 2019. Since the resolution was lower we tried various indexes and bands recommended by Mike Cecil’s research to see what yielded best results. We settled on using the following:

Band 2 (blue)
Band 3 (green)
Band 4 (red)
Band 5 (NIR)
Band 6 (SWIR 1)
Band 7 (SWIR 2)
NDWI index
NWI index
EWI index
Band 5 / Band 4 ratio (substitute for EOMI1)

(chunk 19)

After training the Landsat 8 classifier, we had to redetermine what proportion of the lagoon polygons needed to be classified as lagoon in order to include it in the lagoon count due to the difference in spatial and spectral resolutions between Sentinel-2 and Landsat 8. We tried various proportions to see which produced a lagoon count close to the known number of lagoons in 2019 (235). We settled on classifying any lagoon polygon that had at least 5% of its area classified as a lagoon. At 5%, 239 of the 320 possible lagoon polygons were classified as lagoons in 2019, an overestimation of only 4 lagoons. The bar plot below shows the count of 2019 classified lagoons depending on the proportion cut-off. (chunk 25)

Bar plot comparing the cut-off proportions Used for Landsat 8 2019 classification

The histogram below displays the distribution of lagoon polygon classification for the Landsat 8 classification. The graph demonstrates that most of the lagoon polygons are 0% classified as lagoons. This is in sharp contrast to the distribution of the Sentinel-2 classification displayed after step 10, which classified most of the known lagoons as 100% lagoon. The Landsat 8 classifier therefore misses a lot more lagoons than the Sentinel-2 classifier, which is why we went with a much more lenient cut-off proportion (5% vs 75%).

After settling on a proportion, we created a function that extracts Landsat 8 composites of specified years and returns the number of classified CAFO lagoons for that year. Originally we tried to do it as a for loop to only have to run the code once, but R kept timing out. (chunk 22)

l8_countLagoons() function:

# Function for counting Lagoons with L8
# parameters: 
#   year: year of landsat 8 image
#   proportion: proportion of lagoon area that needs to be classified as lagoon to be counted as lagoon.

# Define function name and parameters
L8_countLagoons <- function(year, proportion) {
  
  ## 1. Import Landsat 8 image for training classifier
  # Landsat 8 data: https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2
  l8_image <- ee$ImageCollection('LANDSAT/LC08/C02/T1_L2')$
    filterBounds(yucatan_feature$geometry())$
    map(function(image){ image$clip(yucatan_feature) })$
    filterDate(paste0(year, '-01-01'), paste0(year, '-12-31'))$
    filter(ee$Filter$lt('CLOUD_COVER', 15))$
    map(maskL8sr)$
    median()$
    setDefaultProjection(crs = "EPSG:32616", scale = 30)
  
  ## 2. Add indexes & ratio to Landsat 8 compound image
  l8_image <- l8_addNWI(l8_image)
  l8_image <- l8_addNDWI(l8_image)
  l8_image <- l8_addEWI(l8_image)
  l8_image <- l8_5_4ratio(l8_image)
  
  ## 3. Classify Landsat 8 images using the trained RandomForest classifier:
  l8_classified <-l8_image$
    select(l8_bands)$                      # use the specified bands for classification
    classify(l8_trainedClassifier)         # use trained clasifier
  
  ## 4. Find average classification inside each lagunas_fc polygon:
  # Add reducer output to the Features in the collection.
  l8_proportionLagoon <- l8_classified$reduceRegions( # only band in the 2019 classified is the presence (0 or 1)
    collection = lagunas_fc,                  # original lagunas fc
    reducer = ee$Reducer$mean(),              # Average
    scale = 30
  )
  
  ## 5. Transfer fc to sf object to get a reliable count.
  l8_proportionLagoon_sf <- ee_as_sf(l8_proportionLagoon)
  
  ## 6. Select Lagoons that are less than 5% classified as lagoon.
  # Select rows with mean <= 0.05 and show only index and mean columns
  l8_lagClass_gte <- l8_proportionLagoon_sf %>%
    filter(mean >= proportion) %>%
    select(index, mean) %>%
    arrange(desc(mean))
  
  ## 7. Find number of lagoons in given year
  l8_num_lagoons <- nrow(l8_lagClass_gte)  # Number of lagoons classified more than 5% as lagoon.
  
  ## 8. Show distribution of % of classified lagoons
  # Create a histogram of the 'mean' column
  l8_lagoonClassDistro <- hist(l8_proportionLagoon_sf$mean, breaks = seq(0, 1, 0.01), 
  xlab = "Avg Classified", ylab = "Frequency", main = paste0("Count of Lagoons included in each classification tier for ", year))
  
  ## 9. Name variables with the year given in argument
  #rm(l8_avglagoon)
  assign(paste0("l8_", year), l8_image, envir = .GlobalEnv)
  assign(paste0("l8_classified_", year), l8_classified, envir = .GlobalEnv)
  assign(paste0("l8_proportionLagoon_sf", year), l8_proportionLagoon_sf, envir = .GlobalEnv)
  assign(paste0("l8_lagClass_gte", proportion, "_sf_", year), l8_lagClass_gte, envir = .GlobalEnv)
  assign(paste0("l8_num_lagoons_", year), l8_num_lagoons, envir = .GlobalEnv)
  assign(paste0("l8_lagoonClassDistro_", year), l8_lagoonClassDistro, envir = .GlobalEnv)
  
  ## 10.Print total lagoon count by year.
  paste0("In ", year, ", there were ", l8_num_lagoons, " CAFO lagoons in the state of Yucatán.")

  ### END
}

We then applied the function to all the years from 2014 to 2022, giving us a lagoon count for each year. We plotted the count of classified lagoons by year to examine possible patterns of CAFO lagoon growth. See bar plot below. (chunk 22, 26)

3. Results (CB)

Lagoon count The CAFO lagoon count for each year is demonstrated below in a table. This table shows the number of lagoons identified in the state of Yucatán. For reference, we originally received a polygon shapefile containing 320 identified lagoons dated from 2000 to 2019.

##           # of CAFOS Identified
## Year_2014                    95
## Year_2015                   124
## Year_2016                   152
## Year_2017                   214
## Year_2018                   185
## Year_2019                   239
## Year_2020                   203
## Year_2021                   166
## Year_2022                   185

l8_countLagoons()

l8_countLagoons() is a function that takes a year from 2014 to present and a percentage. It extracts a Landsat image composite, classifies the image as lagoon or not, makes the classification an sf object, selects the Lagoons that are less than X% classified as lagoon and counts the lagoons classified for that year.

We tested out different proportions of classified lagoons, 5%, 10%, 15%, 20%, 25%, and 33%. For example, if 5% of an existing lagunas_fc polygon was classified as a lagoon in the Landsat 8 image, it would be identified as a lagoon in the classification and counted by the l8_countlagoon() function.

These different proportions of classified lagoons, mean that, for example, in order for a lagoon from the original laguna_fc polygons to be classified as a lagoon in 2014 (or any year), at least 5% of the area must be classified as a lagoon in the 2014 Landsat classification.

The lagoon counts for 2014 - 2022 with a LOESS local regression line showing trend.

Below we compare the first year classification (2014) and lowest lagoon count with the 2019 classification, which had the highest lagoon count.

2014

Map of known 2019 lagoons vs 2014 classified lagoons over the 2014 classification of Yucatán

2019

2019 classifications and lagoons

Color Key:

red = l8-2019 classified lagoons
green = l8-2014 classified lagoons
yellow = S2-2019 classified lagoons
turquoise = Known 2019 lagoons
white = 320 lagoon shapefile
tan = non-lagoons class
blue = lagoons class

4. Discussion (AD)

Our project results have demonstrated, as we expected, a steady increase in CAFO lagoons in Yucatán from 2014 to 2019. To our surprise the number of lagoons dipped in 2020 and 2021. This dip may reflect a decline in CAFO operations due to the COVID-19 global pandemic. However, this may have only been temporary as the count of lagoons increased again in 2022.

R and the rgee package have proven to be a useful tool for automatizing geoprocessing of collections of satellite imagery. Additionally, the ability to access large data files virtually though earth engine servers have allowed for deeper studies of topics’ temporal dimensions.

In regards to our project, there is potential for looking farther back into time with Landsat 7 and especially Landsat 5. This would require the training of new classifiers with older satellite imagery. The lagoon polygon layer would have to be visually re-examined to help verify the accuracy of the classification. This may prove to be time intensive and difficult due to the decreased resolution of older imagery.

Our project unfortunately lacked a validation step due to time constraints. A cross-tabulation of classifications with validation data could help with the fine-tuning of bands, indexes, and ratios used to identify the spectral signatures of the CAFO lagoons. While we were able to almost match the count of known 2019 lagoons with our Landsat 8 classifier, we cannot be confident in the counts generated for other years without better validation of the data.

Finally, our project’s analysis was centered on the researched list of 320 CAFO lagoons. The data provided by Karen Vazques Hudlet, Michael Cecil, and Manuel Llano was essential to our project and we are thankful for their help. It would be interesting to now expand the research beyond the 320 lagoons and try to use remote sensing for identifying new lagoons in order to properly assess the growth of the CAFO industry in the region. As our data showed, CAFO operations may be on the rise again, potentially endangering Yucatán’s ecosystems and communities, there is thus an urgent need for continued attention and research on the topic.

Sources:

FÁBRICAS DE CERDOS EN YUCATÁN. (2022, February 22). ArcGIS StoryMaps. https://storymaps.arcgis.com/collections/3e7203cf44cf417c9b5fe1db7a182293

Mexico - Subnational Administrative Boundaries - Humanitarian Data Exchange. (n.d.). Data.humdata.org. Retrieved May 9, 2023, from https://data.humdata.org/dataset/cod-ab-mex

CAFOs Final Project

André Luiz de Oliveira Domingues and Clio Bate

2023-05-08

Set up code chunks:

Yucatán and CAFO lagoons

1. Introduction (CB)

2. Methods (AD)

3. Results (CB)

4. Discussion (AD)