Matthew Ibarra’s Final Project Was His Favorite

As a SkyTruth Intern, Matthew Ibarra learned new skills and helped protect Native lands.

As I finish up my internship at SkyTruth, I can honestly say that the experience has been everything I imagined it would be and more. My time here was a perfect amalgamation of what I love: namely, an organization that applies technology and gathers and analyzes data to protect the environment. 

When I started my internship at SkyTruth I was unsure of what to expect. I remember the first day I drove into the small town of Shepherdstown, West Virginia. I was worried. For the first time in my life I was working with like-minded individuals with special talents and skills far above my own. I thought that I would have to perform as well as my colleagues right off the bat. However, my fears quickly melted away upon meeting Brendan Jarrell, SkyTruth’s Geospatial Analyst and father to all us interns. Brendan assured me that I would be focusing on my own personal interests, developing practical skills, and applying them to the various projects happening at SkyTruth. Within my first week I became familiar with all the programs I needed for my internship, namely Google Earth Engine and QGIS. Both are programs that are critical in geospatial analysis that were completely new to me, despite having taken Geographic Information System (GIS) courses at West Virginia University. Interning at SkyTruth opened my eyes to the new possibilities of environmental monitoring and I was excited to get started.

My very first day I became familiar with the various satellites that orbit the Earth and provide the imagery that SkyTruth uses on a daily basis. The Landsat and Sentinel satellite missions provide imagery available for free to the public, allowing interns like myself to create maps and interactive data to track activity on Earth. My first task as an intern was to monitor Southeast Asian waters for bilge dumps — oily slicks of wastewater dumped in the ocean by ships. I used Google Earth Engine to access the necessary imagery easily. Then I used QGIS to create the various maps that we post on our Facebook page and blog posts. I found my first bilge dump on February 7, 2020. It was a 35 kilometer slick (almost 22 miles long) off the coast of Malaysia. 

Often, we can identify the likely polluter using Automatic Identification System (AIS) to track vessels at sea. Most vessels constantly send out a radio signal to broadcast their route. When those signal points align with a bilge dump it suggests that the ship is the likely source for that bilge slick. However, not all ships will transmit their signal at all times, and there have even been instances of ships spoofing their signal to different locations. For my first slick I was unable to match a ship’s AIS broadcast to the trail of the bilge dump, but I was able to do so several  times after that. We can’t know for certain who caused this slick, but imagery helps us paint a picture of potential suspects. My first slick pales in comparison to the many slicks I found in later months: later, I captured a few slicks that were over 100 kilometers (more than 62 miles) in length. I was also able to link a ship’s AIS broadcast to the trail of the slick. You can read more about slicks in the South East Asia region in my April 15 blog post here.

 

An example of likely bilge dumping from a vessel identified using Sentinel satellite imagery

Following my introduction to our bilge dumping detection work, I was thrilled to be assigned my first solo project for SkyTruth — updating SkyTruth’s FrackFinder. FrackFinder is an ongoing project at SkyTruth. It aims to keep track of the active oil and natural gas well pads in states such as West Virginia. Drilling permit data is often misleading; sites that are permitted to be drilled may not actually be drilled for several years. In the past, our FrackFinder app was hosted in Carto. Carto is a cloud-based mapping platform that provides limited GIS tools for analysis. I was tasked with giving the application an overhaul and bringing it into Google Earth Engine, a much more powerful and accessible program. 

Learning to code for Earth Engine was challenging for me. I had only one computer science course in college, and that was nearly three years ago. So I was surprised that my first project would revolve around coding. Initially, I was overwhelmed and I struggled to find a place to start. As time went on I slowly became more comfortable with spending large amounts of time solving tiny problems. Brendan was incredibly helpful and patient with teaching me everything I would need to know to be successful. He always made time for me and assisted me with my code numerous times. My finished app is far from perfect but I am proud of the work that I accomplished and I hope that it brings attention to the changing landscape of West Virginia caused by oil and natural gas drilling using hydraulic fracturing (fracking). 

 

The FrackTracker app for West Virginia

My second and final project was creating a visualization about the land surrounding Chaco Culture National Historical Park in New Mexico. Much like the update to the FrackFinder App, it involved the changing landscape surrounding the park due to the increase in fracking. I was tasked with creating a series of still images, an embeddable GIF which shows an animation of the rapid increase in drilling, and an app on Earth Engine that allows the user to zoom in and visually inspect each individual well surrounding the park. In the final months of my internship, I became comfortable using the programs that were foreign to me initially. I created a series of 19 images using QGIS from the years 2000-2018. You can see the collection of images for each year here. SkyTruth’s Geospatial Engineer Christian Thomas assisted me in creating the GIF. 

This project was special to me because I was able to help activists who are advocating for the passage of the Chaco Cultural Heritage Area Protection Act, legislation passed by the U.S. House of Representatives that would effectively create a 10-mile buffer zone surrounding the park and ensure the protection of the area for Native Americans and local communities for generations to come. The Senate has not yet passed the act. When I started my internship at SkyTruth I never would have believed that I would be advocating for protection of Native lands. I always believed issues like these were too big for one person to tackle, but if there’s anything I learned at SkyTruth is that only one person can create real change.

The growth of oil and gas wells within a 15-mile radius of Chaco Culture National Historical Park from 2000 – 2018

After interning at SkyTruth for the past eight months I am happy to say that I feel I have made a difference in the world. I accomplished so much that I thought would be impossible for me initially. I used to think oil slicks were tragedies that happened infrequently, limited to a few times a decade. I was shocked to learn that oily wastewater gets dumped into the ocean so frequently that I was able to capture more than  80 bilge dumps in my eight months at SkyTruth. 

In addition, one of my greatest passions is sustainable energy. I was thrilled to be an advocate for clean energy by showcasing the dangers of an ever-expanding oil and natural gas industry. West Virginia has been my home for the past five years during my time at West Virginia University and I was happy to be able to bring to light one of the growing concerns of the state through the 2018 FrackFinder update. Finally, I was able to advocate for the protection of Native lands through the most meaningful project to me — the Chaco Culture National Historical Park visualizations. It felt incredible fighting for something that was much bigger than myself. As I leave SkyTruth, I will miss contributing to the world in my own way.

SkyTruth has always been more than just a place to intern at for me. I have made unforgettable connections with my colleagues despite the various challenges that we all have to face every day, such as the ongoing COVID-19 pandemic. Never once did I feel that I was alone in my work. I always knew there were people supporting me and encouraging me in my projects even when I was working remotely. I will never forget Christian’s tour of Shepherdstown on my first day or Brendan’s talks about the best Star Wars movie. I cannot thank each of them enough for the patience and kindness they showed me in my short time with them. Everyone at SkyTruth has contributed to my success in some way. I will miss everyone, but I’ll carry my new skills and experiences with me for the rest of my life.   

Drilling Detection with Machine Learning Part 3: Making and mapping predictions

SkyTruth Technical Program Director Ry Covington, PhD explains challenges to generating meaningful predictions from the machine learning model, and outlines solutions.

[This is the final post in a 3-part blog series describing SkyTruth’s effort to automate the detection of oil and gas well pads around the world using machine learning. We hope that this series – and the resources we’ve provided throughout – will help other developers and scientists working on conservation issues to learn from our experience and build their own exciting tools.  You can read the two previous posts in the series here and here.]

SkyTruth Intern Sasha Bylsma and Geospatial Analyst Brendan Jarrell explained how we create training data and implement a machine learning model to detect oil and gas well pads. So, now that we’ve got a trained model, we just need to run it on a few satellite images and put the predictions on a map.  Seems easy enough…  

We started with some Sentinel-2 imagery collected over the Neuquén Basin in central Argentina. This is one of the most heavily drilled areas in Argentina, and we’ve used the Google Earth Engine platform to export a few Sentinel-2 images that we could work with.  

Screenshot of the Neuquén basin in Argentina. 

The images come out of Earth Engine as GeoTIFFs – a pretty standard file format for overhead imagery. We’ve used some Earth Engine magic to reduce the file size of each image so that they’re easier to handle, but they’re still a bit big for the model. The model expects simple, small patches of images: 256 pixels high, 256 pixels wide, and three channels (e.g., Red, Green, Blue) deep. Our Sentinel-2 GeoTIFFs are about 11,000 pixels high by 11,000 pixels wide, so that left us with a few things to figure out:

  • First, the model is expecting small, simple patches of images – no frills, just pixel values. That means that we have to take the geographic information that’s imbedded in the original GeoTIFF and set aside.  So, how do we do that? 
  • Second, how can we evenly slice up the full image into the small patches that the model is expecting?
  • Third, for every small image patch that the model sees, it returns a small, single channel prediction image with values between zero and one. The closer a pixel is to one, the more likely it is that pixel belongs to a drilling site.  But once the model makes predictions on all of the small images, how do we put them back together in the right order to get the original image’s size and shape?
  • Lastly, how do we take the prediction image and convert it into polygons that we can overlay on a map?

These were all things that we’d never done before, so it’s taken us quite a bit of trial and error to figure out how to make things work. In fact, we’re still working on them – we’ve got a workflow in place, but we’re always trying to refine and improve it. For now, let’s just look at what we’ve got working. 

Step 1. Converting our GeoTIFF to a NumPy array

We started off with a pretty small subset of a Sentinel-2 image that we could experiment with. It’s 1,634 pixels high, 1,932 pixels wide, and 3 channels deep. In developer’s parlance, its shape is: (1634, 1932, 3). The image is of Mari Menuco Lake in Argentina. There are a few dozen drilling sites along the southern portion of the lake that seemed like an ideal place to test out our workflow. Once we had everything working like expected, we’d run the larger Sentinel-2 images through.  

First, we used the GDAL Python API to load our image and collect its (a) geotransform and (b) projection. So, what are these two things? Well, basically, the geotransform is the formula that GDAL uses to go from pixel space (think rows and columns) to geographic space (think x [longitude] and y [latitude]), and the projection is just the coordinate reference system of the image. After we had those two pieces of information set aside for later, we pushed all of the image bands into an NumPy array. 

 

# Get geographic information.
projection = tiff.GetProjection()           

# Set spatial reference.
spatial_ref = osr.SpatialReference()     # Create empty spatial reference.
spatial_ref.ImportFromWkt(projection)    # Read the "projection" string.

# Collect all the bands in the .tif image.
bands = [tiff.GetRasterBand(band+1) for band in range(tiff.RasterCount)]

# Read each band as an array.
arrays = [band.ReadAsArray() for band in bands] 

# Combine into a single image. 
image = np.array(arrays)

# Format as (height, width, channels).
image = image.transpose(1,2,0)

GDAL reads and writes images differently than NumPy, so the last thing we did was transpose the axes to put our image in the format that we needed: height, width, channels.   

Step 2. Slicing our NumPy array and running predictions 

The next bit was tricky for us to figure out. We needed to take our image – 1634 by 1932 by 3 (1634, 1932, 3) – and slice it into even squares of (256, 256, 3). Our first problem: neither 1634 nor 1932 divides by 256 evenly, so we needed to figure out how to make the image patches overlap just enough to get a clean division.  

Our second problem: we also needed to keep track of where each patch lived in the larger image so that we could put the predictions back together in the right order later. We ended up giving each patch an ID and collecting the coordinates of their top-left pixel (their minimum x and minimum y). We pushed that information into a pandas dataframe – basically, a 2-D matrix of rows and columns – that we could set aside to rebuild our prediction image later.  

Many thanks to CosmiQ Works and all of the participants in the SpaceNet Challenges; the code snippets and GitHub repos that they’ve made available were immensely helpful for us as we tried to figure out how to implement this step.

 

# Set some variables.
patchSize = 256
overlap = 0.2
height, width, bands = image.shape
imgs, positions = [], []
columns = ['xmin', 'ymin']

# Work through the image and bin it up.
for x in range(0, width-1, int(patchSize*(1-overlap))):    
   for y in range(0, height-1, int(patchSize*(1-overlap))):
      
       # Get top-left pixel.
       xmin = min(x, width-patchSize)
       ymin = min(y, height-patchSize)

       # Get image patch.
       patch = image[ymin:ymin+patchSize, xmin:xmin+patchSize]

       # Set position.
       pos = [xmin, ymin]

       # Add to array.
       imgs.append(patch)
       positions.append(pos)

# Convert to NumPy array.
imageArray = np.array(imgs) / 255.0

# Create position datataframe.
df = pd.DataFrame(positions, columns=columns)
df.index = np.arange(len(positions))

Once we had the new array of patches – 80 patches of 256 by 256 by 3 – it was easy to run the model and generate some predictions.

# And, go. Don't forget the batch dimension.
predictions = model.predict(imageArray, batch_size=20, steps=4)

Step 3. Rebuilding our image

The model returns an array of predictions – (80, 256, 256, 1). The prediction values range from zero to one. So, a pixel value of .82 means that the model is 82% confident that pixel belongs to a drilling site.  

Side by side comparison of an image and its prediction.

We used the pandas dataframe that we made earlier to put all of these predictions back together in the right order and get the original image’s size and shape. The dataframe was where we recorded the ID and top-left pixel (their minimum x and minimum y) of each patch. First, we created an empty image that is the same size and shape as the original. Next, we went through the dataframe, took out each patch, and added it to the new empty image in the right spot (its top-left pixel).  

 

# Create numpy zeros of appropriate shape.
empty_img = np.zeros((height, width, 1))

# Create another zero array to record where pixels get overlaid.
overlays = np.zeros((height, width, 1))

# Iterate through patches.
for index, item in positions.iterrows():

   # Grab values for each row / patch.
   [xmin, ymin] = item

   # Grab the right patch.
   slice = predictions[index]
  
   x0, x1 = xmin, xmin + patchSize
   y0, y1 = ymin, ymin + patchSize

   # Add img_slice to empty_img.
   empty_img[y0:y1, x0:x1] += slice

   # Update count of overlapping pixels.
   overlays[y0:y1, x0:x1] += np.ones((patchSize, patchSize, 1))            

# Normalize the image to get our values between 0 and 1. 
rebuilt_img = np.divide(empty_img, overlay_count)

Most of our team are visual thinkers, so the easiest way for us to imagine rebuilding the image is like covering a blank sheet of white paper in pink sticky-notes, and then smoothing them all down to get a new, blank sheet of pink paper.   

Step 4. Converting our predictions to polygons

After rebuilding our prediction array to be the same size and shape as the original satellite image, we used the GDAL Python API to convert it into polygons that could go on a map. To try and clean things up a bit, we started by selecting only those pixels where the model was more than 70% confident they belonged to a drilling site. We set anything under that threshold to zero. This just helped us to eliminate some noise and clean up the predictions a bit. With that done, we used GDAL to convert the cleaned up prediction image into polygons and reattach the spatial information that we set aside at the beginning (i.e., the geotransform and the projection).  

  

# Band to use.
sourceBand = output.GetRasterBand(1)

# Set up shapefile.
shpDrv = ogr.GetDriverByName("ESRI Shapefile")                                
outFile = shpDrv.CreateDataSource("/gdrive/My Drive/detections.shp")      
layer = outFile.CreateLayer("detections", srs=spatial_ref)                    

# Add field.
idField = ogr.FieldDefn("id", ogr.OFTInteger)               
layer.CreateField(idField)

# And, convert. 
gdal.Polygonize(sourceBand, sourceBand, layer, 0, [], None)

And at this point we had our shapefile. From there, it was easy to upload that shapefile as an asset into our Earth Engine account and have a look at our predictions over satellite imagery. We did a bit of clean up and editing to make sure that all of the polygons look good — what’s sometimes called “human-in-the-loop” modeling — but, for the most part, we were all done.

Screenshot of the polygons in EE.

Lastly, we did a quick assessment to see how well our from-scratch workflow functioned. In the image above, red points are drilling sites that we got correct (true positives), green points are drilling sites that we missed (false negatives), and blue points are places where the model thought there was a drilling site when there wasn’t (false positives). Here are the numbers: 

Total number of validated ground truth points: 239
True Positives: 107
False Positives: 50
False Negatives: 82
Precision: 0.6815286624203821
Recall: 0.5661375661375662
F1-score: 0.6184971098265896

Precision, recall, and F1-score are all just metrics for understanding how a model performs. Your spam folder offers a good example. Imagine that your email’s spam model makes 30 predictions. Precision would measure the percentage of emails flagged as spam that it correctly classified as spam. Recall would measure the percentage of actual spam emails that it correctly classified as spam. Often, these two things are in tension – if the model is more aggressive (i.e., lowers the threshold for an email to be classified as spam) the recall will go up since they’d be capturing more actual spam emails. But the precision would go down, because they’re also capturing more emails overall, and it’s pretty likely that many of those won’t be spam. The F1-score builds off of precision and recall, and it’s probably easiest to think of it as a measure of overall accuracy. In the context of the drilling site work, our precision and recall numbers mean that 68% of the things we’re predicting as drilling sites actually are drilling sites, but we’re only picking up 56% of the drilling sites that are actually out there in the world. 

We hope that this series has helped to demystify the machine learning workflow for others. Figuring out how to build a machine learning data processing pipeline from scratch was an exciting challenge for us. We’re encouraged by our progress so far, but as the metrics above indicate, there’s still lots of room for improvement. We will keep at it because this technology stack — integrating remote sensing with cloud computing and machine learning — is at the heart of our Conservation Vision initiative to automate the analysis of imagery, to solve challenging conservation problems. 

So please stay tuned for updates in the future.  And don’t hesitate to send us a note if you have questions about what we’ve done or ideas about how we can improve things going forward.  We welcome the collaboration!