bellingcat/RS4OSINT
1
0
Fork 0
mirror of https://github.com/bellingcat/RS4OSINT.git synced 2026-06-11 04:58:35 +03:00
Code Issues Packages Projects Releases Wiki Activity

eoghan edits + changed dir structure

Browse Source
This commit is contained in:
Ollie Ballinger
2023-04-17 11:28:29 +01:00
parent 6c5816d659
commit 7ccf85e48e
653 changed files with 44167 additions and 4363 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
chapters/.gitignore vendored Normal file
View File

@@ -0,0 +1 @@
/.quarto/

102
chapters/A2_Remote_Sensing.qmd Normal file
View File

@@ -0,0 +1,102 @@
---
title: Remote Sensing
---
Before learning how to load, process, and analyze satellite imagery in Google Earth Engine, it will be helpful to know a few basic principles of remote sensing. This section provides a brief overview of some important concepts and terminology that will be used throughout the course, including active and passive sensors; spatial, spectral and temporal resolution; and orbits.
## Active and Passive Sensors
[Remote sensing](https://www.sciencedirect.com/topics/medicine-and-dentistry/remote-sensing) is the science of obtaining information about an object or phenomenon without making physical contact with the object. Remote sensing can be done with various types of electromagnetic radiation such as visible, infrared or microwave. The electromagnetic radiation is either emitted or reflected from the object being sensed. The reflected radiation is then collected by a sensor and processed to obtain information about the object.
![](../images/diagram.png)
While most satellite imagery is optical, meaning it captures sunlight reflected by the earth’s surface, Synthetic Aperture Radar (SAR) satellites such as Sentinel-1 work by emitting pulses of radio waves and measuring how much of the signal is reflected back. This is similar to the way a bat uses sonar to “see” in the dark: by emitting calls and listening to echoes.
## Resolution
Resolution is one of the most important attributes of satellite imagery. There are three types of resolution: spatial, spectral, and temporal. Let's look at each of these.
### Spatial Resolution
Spatial resolution governs how "sharp" an image looks. The Google Maps satellite basemap, for example, is really sharp.
Most of the optical imagery that is freely available has relatively low spatial resolution (it looks more grainy than, for example, the Google satellite basemap),
![](../images/Landsat.png)
![](../images/Sentinel2.png)
![](../images/Maxar.png)
### Spectral Resolution
What open access imagery lacks in spatial resolution it often makes up for with *spectral* resolution. Really sharp imagery from MAXAR, for example, mostly collects light in the visible light spectrum, which is what our eyes can see. But there are other parts of the electromagnetic spectrum that we can't see, but which can be very useful for distinguishing between different materials. Many satellites that have a lower spatial resolution than MAXAR, such as Landsat and Sentinel-2, collect data in a wider range of the electromagnetic spectrum.
Different materials reflect light differently. An apple absorbs shorter wavelengths (e.g. blue and green), and reflects longer wavelengths (red). Our eyes use that information-- the color-- to distinguish between different objects. Below is a plot of the spectral profiles of different materials:
<iframe title="Spectral Profiles of Different Materials" aria-label="Interactive line chart" id="datawrapper-chart-b1kcX" src="https://datawrapper.dwcdn.net/b1kcX/3/" scrolling="no" frameborder="0" style="width: 0; min-width: 100% !important; border: none;" height="400"></iframe><script type="text/javascript">!function(){"use strict";window.addEventListener("message",(function(e){if(void 0!==e.data["datawrapper-height"]){var t=document.querySelectorAll("iframe");for(var a in e.data["datawrapper-height"])for(var r=0;r<t.length;r++){if(t[r].contentWindow===e.source)t[r].style.height=e.data["datawrapper-height"][a]+"px"}}}))}();
</script>
The visible portion of the spectrum is highlighted on the left, ranging from 400 nanometers (violet) to 700nm (red). Our eyes (and satellite imagery in the visible light spectrum) can only see this portion of the light spectrum; we can't see UV or infrared wavelengths, for example, though the extent to which different materials reflect or absorb these wavelengths is just as useful for distinguishing between them. The European Space Agency's Sentinel-2 satellite collects spectral information well beyond the visible light spectrum, enabling this sort of analysis. It chops the electromagnetic spectrum up into "bands", and measures how strongly wavelengths in each of those bands is reflected:
![](images/S2_bands.png)
To illustrate why this is important, consider Astroturf (fake plastic grass). Astroturf and real grass will both look green to us, especially from a satellite image. But living plants strongly reflect radiation from the sun in a part of the light spectrum that we can't see (near-infrared). There's a spectral index called the Normalized Difference Vegetation Index (NDVI) which exploits this fact to isolate vegetation in multispectral satellite imagery. So if we look at [Gillette Stadium](https://en.wikipedia.org/wiki/Gillette_Stadium) near Boston, we can tell that the three training fields south of the stadium are real grass (they generate high NDVI values, showing up red), while the pitch in the stadium itself is astroturf (generating low NDVI values, showing up blue).
![VHR image of Gillette Stadium with Sentinel-2 derived NDVI overlay](images/NDVI.jpg)
In other words, even though these fields are all green and indistinguishable to the human eye, their *spectral profiles* beyond the visible light spectrum differ, and we can use this information to distinguish between them.
Astroturf is a trivial example. But suppose we were interested in identifying makeshift oil refineries in Northern Syria that constitute a key source of rents for whichever group controls them. As demonstrated in the ['Refinery Identification'](C2_Refineries.qmd) case study, we can train an algorithm to identify the spectral signatures of oil, and use that to automatically detect them in satellite imagery.
### Temporal Resolution
Finally, the frequency with which we can access new imagery is an important consideration. This is called the **temporal resolution**.
The Google Maps basemap is very high resolution, available globally, and is freely available. But it has no *temporal* dimension: it's a snapshot from one particular point in time. If the thing we're interested in involves *changes* over time, this basemap will be of limited use.
The **"revisit rate"** is the amount of time it takes for the satellite to pass over the same location twice. For example, the Sentinel-2 constellation's two satellites can achieve a revisit rate of 5 days, as shown in this cool video from the European Space Agency:
{{< video https://dlmultimedia.esa.int/download/public/videos/2016/08/004/1608_004_AR_EN.mp4 >}}
Some satellite constellations are able to achieve much higher revisit rates. Sentinel-2 has a revisit rate of 5 days, but SkySat is capable of imaging the same point on earth around 12 times per day! How is that possible? Well, as the video above demonstrated, the Sentinel-2 constellation is composed of two satellites that share the same orbit, 180 degrees apart. In contrast, the SkySat constellation comprises 21 satellites, each with its own orbital path:
{{< video https://assets.planet.com/products/hi-res/Planet_Block_3_HD_1080p.mp4 >}}
This allows SkySat to achieve a revisit rate of 2-3 *hours*. The catch, however, is that you need to pay for it (and it [ain't cheap](https://apollomapping.com/blog/an-update-on-skysat-tasking-pricing-and-video-capabilities)). Below is a comparison of revisit rates for various other optical satellites:
![A chart of revisit times for different satellites from [Sutlieff et. al.(2021)](https://link.springer.com/article/10.1007/s10712-021-09637-5)](images/revisit_chart.png)
## Summary
You should hopefully have a better understanding of what satellite imagery is, and how it can be used to answer questions about the world. In the [next section](A3_Data_acquisition.qmd), we'll look at the various types of satellite imagery stored in the Google Earth Engine catalog.

330
chapters/A3_Data_Acquisition.qmd Normal file
View File

@@ -0,0 +1,330 @@
---
title: Data Acquisition
---
One of the main advantages of GEE is that it hosts several Petabytes of satellite imagery and other spatial data sets, [all in one place](https://developers.google.com/earth-engine/datasets). Among these are many that could prove useful to those investigating illegal mining and logging, estimating conflict-induced damage, monitoring pollution from extractive industries, conducting maritime surveillance without relying on ship transponders, verifying the locations of artillery strikes, tracking missile defense systems and many other topics.
This section highlights ten categories of geospatial data available natively in the GEE catalog, ranging from optical satellite imagery, to atmospheric data, to building footprints. Each sub-section provides an overview of the given data type, suggests potential applications, and lists the corresponding datasets in the GEE catalog. The datasets listed under each heading are **not** an exhaustive list-- there are over 500 in the whole catalog, and the ones listed in this section are simply the ones with the most immediate relevance to open source investigations. If a particular geospatial dataset you want to work with isn't hosted in the GEE catalog, you can upload your own data. We'll cover that in the next section.
## Optical Imagery
![Automatic detection of vehicles using artificial intelligence in high resolution optical imagery. See the [object detection](C5_Object_Detection.qmd) tutorial.](../images/obj_det3.jpg)
Optical satellite imagery is the bread and butter of many open source investigations. It would be tough to list off all of the possible use cases, so here's a handy flowchart:
```{mermaid}
%%{init: {'theme': 'base', 'themeVariables': { 'primaryColor': '#FFFFFF' ,'primaryBorderColor':'#000000' , 'lineColor':'#009933'}}}%%
flowchart
A(Does it happen outside?)
A--> B(Yes)
A--> C(No)
D(Is it very small?)
B-->D
E(Yes)
F(No)
D-->F
D-->E
G(Use optical satellite imagery)
H(Don't use optical satellite imagery)
E-->H
F-->G
C-->H
```
This is, of course, a bit of an exaggeration. But if you're interested in a visible phenomenon that happens outdoors and that isn't very small, chances are an earth-observing satellite has taken a picture of it. What that picture can tell you naturally depends on what you're interested in learning. For a deeper dive into analyzing optical satellite imagery, see the subsection on [multispectral remote sensing.](A2_Remote_Sensing.qmd#multispectral-remote-sensing-remote_sensing).
There are several different types of optical satellite imagery available in the GEE catalog. The main collections are the Landsat and Sentinel series of satellites, which are operated by NASA and the European Space Agency, respectively. Landsat satellites have been in orbit since 1972, and Sentinel satellites have been in orbit since 2015. Norway's International Climate and Forest Initiative (NICFI) has also contributed to the GEE catalog by providing a collection of optical imagery from Planet's PlanetScope satellites. These are higher resolution (4.7 meters per pixel) than Landsat (30m/px) and Sentinel-2 (10m/px), but are only available for the tropics. Even higher resolution imagery (60cm/px) is available from the GEE catalog from the National Agriculture Imagery Program, but it is only available for the United States. For more details, see the "Datasets" section below.
### Applications {.unnumbered}
* Geolocating pictures
- Some of Bellingcat's [earliest work](https://www.bellingcat.com/resources/how-tos/2014/07/09/verification-and-geolocation-tricks-and-tips-with-google-earth/) involved figuring out where a picture was taken by cross-referencing it with optical satellite imagery.
* General surveillance
- [Monitoring](https://web.archive.org/web/20220415054905/https://fas.org/blogs/security/2021/11/a-closer-look-at-chinas-missile-silo-construction/) Chinese missile silo construction.
- Amassing [evidence](https://www.nytimes.com/2022/04/04/world/europe/bucha-ukraine-bodies.html) of genocide in Bucha, Ukraine
* Damage detection
- [Ukraine](https://www.theguardian.com/world/2022/oct/27/before-and-after-satellite-imagery-will-track-ukraine-cultural-damage-un-says)
- [Mali](https://reliefweb.int/report/mali/satellite-imagery-conflict-affected-areas-how-technology-can-support-wfp-emergency)
- [Around the World](https://www.pnas.org/doi/pdf/10.1073/pnas.2025400118)
* Verifying the locations of artillery/missile/drone strikes
- The [2019 attack](https://www.cnbc.com/2019/09/17/satellite-photos-show-extent-of-damage-to-saudi-aramco-plants.html) on Saudi Arabia's Abqaiq oil processing facility.
* Monitoring illegal mining/logging
- Global Witness [investigation](https://www.globalwitness.org/en/campaigns/natural-resource-governance/myanmars-poisoned-mountains/) into illegal mining by militias in Myanmar.
- Tracking [illegal logging](https://www.theguardian.com/environment/2016/mar/02/new-satellite-mapping-a-game-changer-against-illegal-logging) across the world.
### Datasets {.unnumbered}
| Sensor | Timeframe | Resolution | Coverage |
| ----------- | ------------ | ---------- | -------- |
| [Landsat 1-5](https://developers.google.com/earth-engine/datasets/catalog/landsat-mss) | 1972–1999 | 30m | Global |
| [Landsat 7](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LE07_C02_T1_L2) | 1999–2021 | 30m | Global |
| [Landsat 8](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2) | 2013–Present | 30m | Global |
| [Landsat 9](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC09_C02_T1_L2) | 2021–Present | 30m | Global |
| [Sentinel-2](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR_HARMONIZED) | 2015–Present | 10m | Global |
| [NICFI](https://developers.google.com/earth-engine/datasets/tags/nicfi) | 2015-Present | 4.7m | Tropics |
| [NAIP](https://developers.google.com/earth-engine/datasets/catalog/USDA_NAIP_DOQQ) | 2002-2021 | 0.6m | USA |
## Radar Imagery
![Ships and interference from a radar system are visible in Zhuanghe Wan, near North Korea.](../images/radar%20ships.jpg)
Synthetic Aperture Radar imagery (SAR) is a type of remote sensing that uses radio waves to detect objects on the ground. SAR imagery is useful for detecting objects that are small, or that are obscured by clouds or other weather phenomena. SAR imagery is also useful for detecting objects that are moving, such as ships or cars.
### Applications {.unnumbered}
* Change/Damage detection
* Tracking military radar systems
* Maritime surveillance
* Monitoring illegal mining/logging
### Datasets {.unnumbered}
| Sensor | Timeframe | Resolution | Coverage |
| ----------- | ------------ | ---------- | -------- |
| [Sentinel 1](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD) | 2014-Present | 10m | Global |
## Nighttime Lights
![A timelapse of nighttime lights over Northern Iraq showing the capture and liberation of Mosul by ISIS.](../images/Figure_1.gif)
Satellite images of the Earth at night are a useful proxy for human activity. The brightness of a given area at night is a function of the number of people living there and the nature of their activities. The effects of conflict, natural disasters, and economic development can all be inferred from changes in nighttime lights.
The timelapse above reveals a number of interesting things: The capture of Mosul by ISIS in 2014 and the destruction of its infrastructure during the fighting (shown as the city darkening), as well as the liberation of the city by the Iraqi military in 2017 are all visible in nighttime lights. The code to create this gif, as well as a more in-depth tutorial on the uses of nighttime lights, can be found in the ["War at Night"](C1_Lights.qmd) case study.
### Applications {.unnumbered}
* Damage detection
* Identifying gas flaring/oil production
* Identifying urban areas/military bases illuminated at night
### Datasets {.unnumbered}
| Sensor | Timeframe | Resolution | Coverage |
| ----------- | ------------ | ---------- | -------- |
| [DMSP-OLS](https://developers.google.com/earth-engine/datasets/catalog/NOAA_DMSP-OLS_NIGHTTIME_LIGHTS) | 1992-2014 | 927m | Global |
| [VIIRS](https://developers.google.com/earth-engine/datasets/catalog/NOAA_VIIRS_DNB_MONTHLY_V1_VCMSLCFG) | 2014-Present | 463m | Global |
## Climate and Atmospheric Data
![Sulphur Dioxide plume resulting from ISIS attack on the Al-Mishraq Sulphur Plant in Iraq](../images/mishraq_small.gif){width=100%}
Climate and atmospheric data can be used to track the effects of conflict on the environment. The European Space Agency's Sentinel-5p satellites measure the concentration of a number of atmospheric gasses, including nitrogen dioxide, methane and ozone. Measurements are available on a daily basis at a fairly high resolution (1km), allowing for the detection of localized sources of pollution such as oil refineries or power plants. For example, see this [Bellingcat article](https://www.bellingcat.com/resources/2021/04/15/what-oil-satellite-technology-and-iraq-can-tell-us-about-pollution/) in which Wim Zwijnenburg and I trace pollution to specific facilities operated by multinational oil companies in Iraq.
The Copernicus Atmosphere Monitoring Service (CAMS) provides similar data at a lower spatial resolution (45km), but measurements are available on an hourly basis. The timelapse above utilizes CAMS data to show a sulfur dioxide plume resulting from an ISIS attack on the Al-Mishraq Sulphur Plant in Iraq. The plant was used to produce sulphuric acid, for use in fertilizers and pesticides. The attack destroyed the plant, causing a fire which burned for a month and released [21 kilotons](https://earthobservatory.nasa.gov/images/88994/sulfur-dioxide-spreads-over-iraq) of sulfur dioxide into the atmosphere per day; the largest human-made release of sulfur dioxide in history.
### Applications {.unnumbered}
* Monitoring of airborne pollution
* Tracing pollution back to specific facilities and companies
* Visualizing the effects of one-off environmental catastrophes
- Nordstream 1 leak
- ISIS setting Mishraq sulphur plant on fire
### Datasets {.unnumbered}
| Sensor | Timeframe | Resolution | Coverage |
| ----------- | ------------ | ---------- | -------- |
| [CAMS NRT](https://developers.google.com/earth-engine/datasets/catalog/ECMWF_CAMS_NRT) | 2016-Present | 44528m | Global |
| [Sentinel-5p](https://developers.google.com/earth-engine/datasets/catalog/sentinel-5p) | 2018-Present | 1113m | Global |
## Mineral Deposits
![Zinc deposits across Central Africa](../images/mining.jpg)
Mining activities often play an important role in conflict. According to an influential [study](https://www.aeaweb.org/articles?id=10.1257/aer.20150774), "the historical rise in mineral prices might explain up to one-fourth of the average level of violence across African countries" between 1997 and 2010. Data on the location of mineral deposits can be used to identify areas where mining activities are likely to be taking place, and several such datasets are available in Google Earth Engine.
### Applications {.unnumbered}
* Monitoring mining activity
* Identifying areas where mining activities are likely to be taking place
* Mapping the distribution of resources in rebel held areas in conflicts fueled by resource extraction
### Datasets {.unnumbered}
| Sensor | Timeframe | Resolution | Coverage |
| ----------- | ------------ | ---------- | -------- |
| [iSDA](https://developers.google.com/earth-engine/datasets/tags/isda) | 2001-2017 | 30m | Africa |
## Fires
![Detected fires over Ukraine since 27/02/2022 showing the frontline of the war](../images/fires.jpg)
Earth-observing satellites can detect "thermal anomalies" (fires) from space. NASA's Fire Information for Resource Management System (FIRMS) provides daily data on active fires in near real time, going back to the year 2000. Carlos Gonzales wrote a comprehensive [Bellingcat article](https://www.bellingcat.com/resources/2022/10/04/scorched-earth-using-nasa-fire-data-to-monitor-war-zones/) on the use of FIRMS to monitor war zones from Ukraine to Ethiopia. The map above shows that FIRMS detected fires over Eastern Ukraine trace the frontline of the war.
FIRMS data are derived from the MODIS satellite, but only show the central location and intensity of a detected fire. Another MODIS product (linked in the table below) generates a monthly map of burned areas, which can be used to assess the spatial extent of fires.
### Applications {.unnumbered}
* Identification of possible artillery strikes/fighting in places like Ukraine
* Environmental warfare and "scorched earth" policies
* Large scale arson
- e.g. [Refugee camps burned down in Myanmar](https://citizenevidence.org/2021/02/26/using-viirs-fire-data-for-human-rights-research/)
### Datasets {.unnumbered}
| Sensor | Timeframe | Resolution | Coverage |
| ----------- | ------------ | ---------- | -------- |
| [FIRMS](https://developers.google.com/earth-engine/datasets/catalog/FIRMS) |2000-Present | 1000m | Global |
| [MODIS Burned Area](https://developers.google.com/earth-engine/datasets/catalog/CIESIN_GPWv411_GPW_Population_Count) | 2000-Present | 500m | Global |
## Population Density Estimates
![Population density estimates around Pyongyang, North Korea](../images/pop.jpg)
Sometimes, we may want to get an estimate of the population in a specific area to ballpark how many people might be affected by a natural disaster, a counteroffensive or a missile strike. You can't really Google "what is the population in this rectangle I've drawn in Northeastern Syria?" and get a good answer. Luckily, there are several spatial population datasets hosted in GEE that let you do just that. Some, such as WorldPop, provide estimated breakdowns by age and sex as well. However, it is extremely important to bear in mind that these are **estimates**, and will **not** take into account things like conflict-induced displacement. For example, Oak Ridge National Laboratory's LandScan program has released high-resolution population data for Ukraine, but this pertains to the pre-war population distribution. The war has radically changed this distribution, so these estimates no longer reflect where people *are*. Still, this dataset could be used to roughly estimate displacement or the number of people who will need new housing.
### Applications: {.unnumbered}
* Rough estimates of civilians at risk from conflict or disaster, provided at a high spatial resolution
### Datasets {.unnumbered}
| Sensor | Timeframe | Resolution | Coverage |
| ----------- | ------------ | ---------- | -------- |
| [Worldpop](https://developers.google.com/earth-engine/datasets/tags/worldpop) |2000-2021 | 92m | Global |
| [GPW](https://developers.google.com/earth-engine/datasets/catalog/CIESIN_GPWv411_GPW_Population_Count) | 2000-2021 | 927m | Global |
| [LandScan](https://developers.google.com/earth-engine/datasets/catalog/DOE_ORNL_LandScan_HD_Ukraine_202201) | 2013–Present | 100m | Ukraine |
## Building Footprints
![Building footprints in Mariupol, Ukraine colored by whether the building is damaged](../images/footprints.png)
A building footprint dataset contains the two dimensional outlines of buildings in a given area. Currently, GEE hosts one building footprint dataset which covers all of Africa. In 2022, Microsoft released a free [global building footprint dataset](https://www.microsoft.com/en-us/maps/building-footprints), though to use it in Earth Engine you'll have to download it from their [GitHub page](https://github.com/Microsoft/USBuildingFootprints) and upload it manually to GEE. The same goes for OpenStreetMap (OSM), a public database of building footprints, roads, and other features that also contains useful annotations for many buildings indicating their use. [Benjamin Strick](https://www.youtube.com/watch?v=bJkV3l5Haq0) has a great youtube video on conducting investigations using OSM data.
### Applications: {.unnumbered}
* Joining damage estimate data with the number of buildings in an area
### Datasets {.unnumbered}
| Dataset | Timeframe | Coverage |
| ----------- | ------------ | -------- |
| [Open Buildings](https://developers.google.com/earth-engine/datasets/catalog/GOOGLE_Research_open-buildings_v2_polygons) |2022 | Africa |
## Administrative Boundaries
![Second-level administrative boundaries in Yemen](../images/fao_gaul.jpg)
Spatial analysis often has to aggregate information over a defined area; we may want to assess the total burned area by province in Ukraine, or count the number of Saudi airstrikes by district in Yemen. For that, we need data on these administrative boundaries. GEE hosts several such datasets at the country, province, and district (or equivalent) level.
### Applications {.unnumbered}
* Quick spatial calculations for different provinces/districts in a country
- e.g. counts of conflict events by district over time
### Datasets {.unnumbered}
| Dataset | Timeframe | Coverage |
| ----------- | ------------ | -------- |
| [FAO GAUL](https://developers.google.com/earth-engine/datasets/tags/gaul) |2015 | Global |
## Global Power Plant Database
![Power plants in Ukraine colored by type](../images/power.jpg)
The Global Power Plant Database is a comprehensive, open source database of power plants around the world. It centralizes power plant data to make it easier to navigate, compare and draw insights. Each power plant is geolocated and entries contain information on plant capacity, generation, ownership, and fuel type. As of June 2018, the database includes around 28,500 power plants from 164 countries. The database is curated by the [World Resources Institute (WRI)](https://datasets.wri.org/dataset/globalpowerplantdatabase).
### Applications: {.unnumbered}
* Analyzing the impact of conflict on critical infrastructure.
- e.g. fighting in Ukraine taking place around nuclear power facilities.
* Could be combined with the atmospheric measurements of different pollutants and the population estimates data to assess the impact of various forms of energy generation on air quality and public health.
### Datasets {.unnumbered}
| Dataset | Timeframe | Coverage |
| ----------- | ------------ | -------- |
| [GPPD](https://developers.google.com/earth-engine/datasets/catalog/WRI_GPPD_power_plants) |2018 | Global |

1576
chapters/B1_Getting_Started.qmd Normal file
View File

File diff suppressed because it is too large Load Diff

1320
chapters/B2_Interpreting_Images.qmd Normal file
View File

File diff suppressed because it is too large Load Diff

3121
chapters/B3_Image_Series.qmd Normal file
View File

File diff suppressed because it is too large Load Diff

1949
chapters/B4_Vectors_Tables.qmd Normal file
View File

File diff suppressed because it is too large Load Diff

296
chapters/C1_Lights.qmd Normal file
View File

@@ -0,0 +1,296 @@
---
title: War at Night
format:
html:
number-sections: true
number-offset: [1,5]
---
Satellite images of Syria taken at night capture a subtle trace left by human civilization: lights. Apartment buildings, street lights, highways, power plants -- all are illuminated at night and can be seen from space. Researchers often use these nighttime lights signatures to track development; as cities grow, villages receive power and infrastructure is built, areas emit more light. But this works both ways. As cities are demolished, villages burned and highways cutoff, they stop emitting lights.
In this tutorial, we'll use satellite images of Iraq taken at night to track the destruction caused by the fight against the Islamic State. We'll use the VIIRS nighttime lights dataset, which is a collection of satellite images taken by the Visible Infrared Imaging Radiometer Suite (VIIRS) on the Suomi NPP satellite. VIIRS is a sensor that can detect light in the visible and infrared spectrum, and is capable of taking images at night. A link to the GEE code for this section can be found [here](https://code.earthengine.google.com/2cf77d8cb9afd76b73100637fbffdf5d).
## Pre-Processing
First, let's start by importing a few useful packages written by [Gennadii Donchyts](https://twitter.com/gena_d). We'll use `utils` and `text` to annotate the date of each image on the timelapse. We'll also define an Area of Interest (AOI), which is just a rectangle. You can do this manually by clicking the drawing tools in the top left. I've drawn an AOI over the area covering Mosul, Irbil and Kirkuk in Northern Iraq.
``` js
var utils = require("users/gena/packages:utils");
var text = require("users/gena/packages:text");
// define the Area of Interest (AOI)
var AOI = ee.Geometry.Polygon(
[[[42.555362833405326, 36.62010778397765],
[42.555362833405326, 35.18296243288332],
[44.681217325592826, 35.18296243288332],
[44.681217325592826, 36.62010778397765]]])
// start and end dates for our gif
var startDate = '2013-01-01';
var endDate = '2018-01-01';
// a filename for when we export the gif
var export_name='qayyarah_viirs'
// A palette to visualize the VIIRS imagery. This one is similar to Matplotlib's "Magma" palette.
var viirs_palette = [
"#000004",
"#320a5a",
"#781b6c",
"#bb3654",
"#ec6824",
"#fbb41a",
"#fcffa4",
];
// Visualisation parameters for the VIIRS imagery, defining a minimum and maximum value, and referencing the palette we just created
var VIIRSvis = { min: -0.1, max: 1.6, palette: viirs_palette };
```
Next, we'll load the VIIRS nighttime lights imagery. We want to select the `avg_rad` band of the image collection, and filter blank images. Sometimes, we get blank images over an area in VIIRS if our AOI is on the edge of the satellite's imaging swath. We can filter these images, similarly to how we filter for cloudy images in Sentinel-2:
``` js
var VIIRS= ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG")
.select('avg_rad')
// Calculate the sum of the 'avg_rad' band within the AOI
.map(function(image) {
var blank=image.reduceRegions({ // reduceRegions is a function that allows us to reduce the values of a band within a
collection: AOI, // geometry. In this case, we're reducing the values of the 'avg_rad' band within the AOI
reducer: ee.Reducer.sum(), // We're using the sum reducer, which will sum the values of the 'avg_rad' band
scale: 10}) // We're reducing the values of the 'avg_rad' band at a scale of 10m
.first() // We only want the first element of the collection, which is the sum of the 'avg_rad' band within the AOI
.get('sum') // We want the value of the 'sum' property, which is the sum of the 'avg_rad' band within the AOI
// For each image, define a property 'blank' that stores the sum of the 'avg_rad' band within the AOI.
// We're also going to take a base 10 log of the image-- this will help us visualize the data by dampening extreme values
return image.set('blank', blank).log10().unmask(0)
})
// Now, we can filter images which are fully or partially blank over our AOI
.filter(ee.Filter.gt('blank', 10))
// Finally, we filter the collection to the specified date range
.filterDate(startDate, endDate)
```
Let's have a look at the first image in the collection to make sure everything's looking right. We'll set the basemap to satellite and center our AOI:
``` js
Map.setOptions('HYBRID')
Map.centerObject(AOI)
Map.addLayer(VIIRS.first(),VIIRSvis,'Nighttime Lights')
```
![](../images/iraq_check.png)
If we decrease the opacity of the VIIRS layer, we can see the cities of Mosul, Erbil and Kirkuk shining brightly at night. We can also see a string of bright lights between Kirkuk and Erbil -- these are methane flares from oil wells.
## Analysis
Having pre-processed the VIIRS imagery, we can now define a function `gif` that will take:
1. An image collection (`col`, in this case the nighttime lights imagery `VIIRS`)
2. Visualization parameters (`col_vis`, in this case `VIIRSvis`)
3. An Area of Interest `AOI`
The function will then return a timelapse.
``` js
var gif = function (col, col_vis, AOI) {
// Define the date annotations to be printed in the top left of the gif in white
var annotations = [
{
textColor: "white",
position: "left",
offset: "1%",
margin: "1%",
property: "label",
// Dynamically size the annotations according to the size of the AOI
scale: AOI.area(100).sqrt().divide(200),
},
];
// Next, we want to map over the image collection,
var rgbVis = col.map(function (image) {
// Get the date of the image and format it
var start = ee.Date(image.get("system:time_start"));
var label = start.format("YYYY-MM-dd");
// And visualize the image using the visualization parameters defined earlier.
// We also want to set a property called "label" that stores the formatted date
return image.visualize(col_vis).set({ label: label });
});
// Now we use the label property and the annotateImage function from @gena_d to annotate each image with the date.
rgbVis = rgbVis.map(function (image) {
return text.annotateImage(image, {}, AOI, annotations);
});
// Define GIF visualization parameters.
var gifParams = {
maxPixels: 27017280,
region: AOI,
crs: "EPSG:3857",
dimensions: 640,
framesPerSecond: 5,
};
// Export the gif to Google Drive
Export.video.toDrive({
collection: rgbVis, // the image collection
description: export_name, // the name of the file
dimensions: 1080, // the dimensions of the gif
framesPerSecond: 5, // the number of frames per second
region: AOI, // the area of interest
});
// Print the GIF URL to the console.
print(rgbVis.getVideoThumbURL(gifParams));
// Render the GIF animation in the console.
print(ui.Thumbnail(rgbVis, gifParams));
};
```
Ok that was a pretty big chunk of code. But the good news is that we basically never have to touch it again, since we can just feed it different inputs. For example, if I want to generate a gif of night time lights over a different area, it's as simple as dragging the AOI. If I want to look at a different time period, I can just edit the `startDate` and `endDate` variables. And if I want to visualize an entirely different type of satellite imagery -- Sentinel-1, Sentinel-2, or anything else, all I have to do is change the image collection (`col`) and visualization parameters (`col_vis`) variables. Now, let's look at some timelapses.
### The Fall of Mosul
The function returns a timelapse of nighttime lights over Northern Iraq:
``` js
gif(VIIRS, VIIRSvis, AOI);
```
![I've done a bit of post-processing to this gif, adding more annotations and blending between frames to make it a bit smoother. I typically use [ffmpeg](https://ffmpeg.org/) and [ezgif](https://ezgif.com/) for the finishing touches. ](../images/Figure_1.gif)
This timelapse gives a play-by-play of one of the most important campaigns in the war against the Islamic State. In the first few frames, Mosul is under the control of the Kurdistan Regional Government (KRG). In the summer of 2014, ISIS captures the city, and power is cut off. Mosul and many villages along the Tigris river are plunged into darkness. In 2015, the front line in the campaign to retake the city emerges around Mosul, advancing in 2016 and 2017. Mosul is eventually retaken by the KRG in 2017, after which it brightens once again as electricity is restored.
### The Qayyarah Fires
Farther south, there is an interesting detail. Above the "h" in "Qayyarah", a bright set of lights emerges just before Mosul is recaptured, around December 2016. Fleeing Islamic State fighters [set fire to the Qayyarah oilfields](https://time.com/iraq-fires/), which burned for months.
Using the VIIRS data we've already loaded, we can further analyze the effect of the conflict using a chart. First, let's define two rectangles (again, you can draw these) over Mosul and Qayyarah:
``` js
var mosul = ee.Feature(
ee.Geometry.Polygon(
[[[43.054977780266675, 36.438274276521234],
[43.054977780266675, 36.290642221212416],
[43.24792516796199, 36.290642221212416],
[43.24792516796199, 36.438274276521234]]], null, false),
{
"label": "Mosul",
"system:index": "0"
}),
qayyarah = ee.Feature(
ee.Geometry.Polygon(
[[[43.08240275545117, 35.8925587996721],
[43.08240275545117, 35.77899970860588],
[43.26642375154492, 35.77899970860588],
[43.26642375154492, 35.8925587996721]]], null, false),
{
"label": "Qayyarah",
"system:index": "0"
})
// Let's put these together in a list
var regions=[qayyarah, mosul]
```
Once we've got the rectangles, we can make a chart that will take the mean value of the VIIRS images in each rectangle over time:
``` js
var chart =
ui.Chart.image
.seriesByRegion({
imageCollection: VIIRS,
regions: regions,
reducer: ee.Reducer.mean(),
seriesProperty:'label'
}).setOptions({
title: 'Nighttime Lights'
});
print(chart)
```
![](../images/qayyarah_chart.png)
We can clearly see Mosul (the red line) darkening in 2014 as the city is taken by ISIS. During this period the Qayyarah oil fields are, as we might expect, quite dark. All of a sudden in 2016 Qayyarah becomes brighter at night than the city of Mosul ever was, as the oilfields are set on fire. Then, almost exactly when the blaze in Qayyarah is extinguished and the area darkens (i.e. when the blue line falls back to near zero), Mosul brightens once again (i.e. the red line rises) as the city is liberated.
<!--
### The Battle for Aleppo
The images below were taken between 2012 and 2014. Vast swaths of the city darken as neighborhoods are razed by fighting.
<timelapse>
Though this is a trend that can be observed across the country, nowhere is the decline in nightlights more visible than in Aleppo. Below is a comparison of longitudinal trends in nightlights signatures between several cities:
<graph>
The most salient trend is Aleppo plummeting over the course of 2012, and becoming steadily darker over the course of the next four years. Raqqa drops in 2012 as well, but remains in flux until 2017, when the battle to reclaim the city plunges it into near total darkness. Damascus also experiences a dip in 2012, but stabilizes relatively quickly. The Turkish city of Gaziantep -- less than 100km from Aleppo and roughly 1/5th the size -- stands in stark contrast to the Syrian cities, becoming progressively brighter over the entire period.
Another interesting pattern here is the difference in seasonal trends in nightlights. Under normal circumstances in this part of the world, cities become brighter at night during the summer months. Restaurants, bars, and markets stay open later and conduct business outdoors. Gaziantep, which still attracts scores of tourists every year, displays pronounced seasonality. Damascus, the most stable of the three Syrian cities, also maintains a seasonal trend throughout the war. In contrast, both Raqqa and Aleppo maintain extremely low and roughly constant levels of nightlights year-round during the periods following intense fighting.
Reliable economic data for Syria haven't been available for nearly a decade, and assessing the country's recovery is consequently difficult. But subtle indications of economic growth are visible above: all three Syrian cities have been on a steady upward trend since 2017, and beginning to display seasonal variation once again. -->
<!-- ### Fighting for Oil
Throughout the war, sudden massive spikes in nightlights signatures can be observed throughout the country. In the center of the map just west of Palmyra, some particularly large spikes occur in 2017:
These flashes of light show gas wells being set on fire, a common form of sabotage carried out by retreating Islamic State fighters. Modified Sentinel-2 imagery of the Hayyan gas field (indicated by the green box above) shows this in greater detail. Substituting the Red band in an RGB image with Near Infrared (NIR) highlights thermal signatures, showing fires burning brightly even during the day.
The large complex on the right is the Hayyan Gas Plant, which produced nearly one third of Syria's electricity. The plant and its associated wells changed hands several times throughout the war, but were under Islamic State control until February 2017. In the video below, Islamic State fighters can be seen rigging the plant with explosives and destroying it on January 8th:
In February, three Russian oil and gas companies (Zarubij Naft, Lukoil and Gazprom Neft) were given restoration, exploration and production rights to the hydrocarbon deposits West of Palmyra. On January 12th, 2017, the Syrian Army's 5th Legion and Russian special forces launched a counterattack known as the "Palmyra offensive", with the aim of retaking several important hydrocarbon deposits including Hayyan.
The timing of well fires aligns closely with a detailed timeline of the campaign.The Near Infrared Sentinel-2 image below shows the layout of the Hayyan Gas Plant and the wells in the Hayyan gas field:
The Syrian Army took the Hayyan gas field on [February 4th](https://www.almasdarnews.com/article/syrian-army-liberates-hayyan-gas-fields-west-palmyra/), and retreating ISIS fighters set fire to wells 1, and 3. However, ISIS managed to briefly retake the Hayyan field on [February 7th](https://www.almasdarnews.com/article/isis-retakes-hayyan-gas-fields-new-bid-expand-west-palmyra/), setting fire to wells 2 and 4. These moments in the Palmyra Offensive are captured in NIR signatures
Interestingly, despite the massive explosion caused by the bombing of the Hayyan Gas Plant, no prolonged thermal anomalies were detected over the area of the plant itself. The well fires, on the other hand, lasted for months. Below is an image of well fire at the Hayyan field taken from this [video](https://www.youtube.com/watch?v=WFe9abYyqK0); based on the nearby infrastructure and date (04/02/2017) of posting, it is likely Well-3.
-->

355
chapters/C2_Refineries.qmd Normal file
View File

@@ -0,0 +1,355 @@
---
title: Refinery Identification
---
*Topics: multispectral satellite imagery, machine learning, informal economies, war.*
In Syria, over a decade of war has ravaged one of its most important industries. Oil is a basic necessity for local residents who need to heat their homes and keep the lights on. It's also an important source of income for armed groups who control production; by some estimates the Islamic State was making over [$40 million](https://www.rand.org/blog/2017/10/oil-extortion-still-paying-off-for-isis.html) per month in oil revenues, and the Syrian Democratic Forces were making [$120 million](https://www.al-monitor.com/originals/2021/08/syrian-government-kurds-discuss-plans-oil-trade) per year selling oil to their adversaries, the Syrian Government. This, in turn, has made oil production facilities a frequent target of [airstrikes](https://www.gov.uk/government/publications/british-forces-air-strikes-in-iraq-monthly-list/january-2015), leading to catastrophic environmental consequences.
The destruction of Syria's oil infrastructure and its importance as a source of revenue for armed groups has led to a massive rise in makeshift oil extraction and refining. These makeshift refineries are often constructed by digging a large pit, lining it with a tarp and filling it with polluted water. A furnace heats crude oil, which is run through a pipe cooled by the basin and collected in drums:
![credit: PAX/Wim Zwijnenburg](images/makeshift-refining.png)
Wim Zwijnenburg wrote an excellent [Bellingcat article](https://www.bellingcat.com/news/2020/04/24/dying-to-keep-warm-oil-trade-and-makeshift-refining-in-north-west-syria/) on the subject, which you should read before going any further in this tutorial. In the article, Wim notes that these facilities "can be spotted by the ditch and the black spot at the end with oil waste residues, which blacken the soil around the furnace." These refineries also frequently leak, blackening large swaths of land around them.
![source: https://www.bellingcat.com/news/2020/04/24/dying-to-keep-warm-oil-trade-and-makeshift-refining-in-north-west-syria/](images/refinery.png)
Looking around Northwestern Syria, we can see agricultural fields pockmarked by these makeshift refineries (you can pan around and zoom in):
```{python}
#| echo: false
from ipyleaflet import Map, Marker, basemaps, basemap_to_tiles
m = Map(
basemap=basemap_to_tiles(
basemaps.Esri.WorldImagery
),
center=(36.936622, 42.118185),
zoom=14
)
m
```
Previous efforts to quantify informal oil production have involved manually sifting through satellite imagery and counting them. This is a painstaking process that leaves a number of important questions unanswered. Even if we were to count all of the individual refineries, could we get an estimate of the polluted area? What if we wanted to count the refineries in a new part of Syria? Or get annual or even monthly estimates of new refineries?
Below is an Earth Engine application that automates the detection of makeshift refineries in Northeastern Syria, using multispectral satellite imagery and machine learning. Blue dots represent the locations of predicted makeshift oil refineries and general oil pollution, while red areas indicate areas predicted to be contaminated by oil.
:::{.column-page}
<iframe src='https://ollielballinger.users.earthengine.app/view/rojavaoil' width='100%' height='700px'></iframe>
:::
You can draw an Area of Interest (AOI) and get the total number of contaminated points as well as the total number of contaminated square meters within the AOI. Drawing multiple AOIs will show a running total of these statistics. It's not perfect -- it misses some refineries and falsely identifies some others -- but it is generally quite accurate; you can visually verify the results of the prediction by zooming in using the "+" button. You can toggle different layers using the "layers" tab as well. This tool could be used to get an estimate of oil production in a user-defined area, and eventually to direct cleanup efforts. The fullscreen version of the application can be found [here](https://ollielballinger.users.earthengine.app/view/rojavaoil), and the source code [here](https://code.earthengine.google.com/7a80f10412e1eb2a4d2c5d95989e70bd). This tutorial will first cover the basics of multispectral remote sensing, before moving into a step-by-step guide in the construction of this model.
# Machine Learning Workflow
## Pre-Processing
As always, the first step in our project will be to load and pre-process satellite imagery. For this project, we'll be using Sentinel-2 imagery. Let's load imagery from 2020-2021, filter out cloudy images, and define visualization parameters:
```js
var start='2020-04-01'
var end='2021-07-01'
var bands = ['B2', 'B3', 'B4','B5','B6','B7','B8', 'B8A','B11','B12']
var sentinel = ee.ImageCollection('COPERNICUS/S2_SR')
.filter(ee.Filter.date(start, end))
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))
.mean()
.select(bands)
var s_rgb = {
min: 0.0,
max: 3000,
bands:['B4', 'B3', 'B2'],
opacity:1
};
```
When loading the Sentinel-2 imagery, I've also only selected the bands that we will ultimately use in our analysis. There are a number of other bands included in the data that we don't need. I've omitted a few bands (B1, B9, B10) because they're collected at a much lower spatial resolution (60 meters) compared to the other bands.
A couple types of landcover are so readily identifiable that we can remove them with thresholds. Water and vegetation both have spectral indices; we looked at NDVI above, but there's a similar one for water called NDWI. These can be calculated from Sentinel-2 imagery as follows:
```js
var ndvi=sentinel.normalizedDifference(['B8','B4'])
.select(['nd'],['ndvi'])
var ndwi=sentinel.normalizedDifference(['B3','B8'])
.select(['nd'],['ndwi'])
```
We use the `normalizedDifference` function and specify which bands we want to use for each index. NDVI uses the red and near infrared bands (B4 and B8), while NDWI uses bands 3 and 8. Finally, we want to rename the resulting band from 'nd' to the name of the spectral index.
Now we can use these indices to filter out water and vegetation. We do this using the `updateMask` function, and specify that we want to remove areas that have an NDVI value higher than 0.2 and anNDWI value higher than 0.3. You can play around with these thresholds until you achieve the desired results.
```js
var image=sentinel.updateMask(ndwi.lt(0.3))
.updateMask(ndvi.lt(0.2))
.addBands(ndvi)
.select(bands)
```
We also want to only select bands that are relevant to our analysis; Sentinel
Finally, let's clip the image to our Area of Interest (AOI) and add it to the map using the visualization parameters we defined earlier.
```js
Map.addLayer(image.clip(AOI), s_rgb, 'Sentinel');
```
![water and vegetation have been removed from this Sentinel-2 image. What remains is largely fallow agricultural land, urban areas, and oil spills.](images/rojava_preprocessed.png)
Now that we've loaded and pre-porcessed our satellite imagery, we can proceed with the rest of our task. Ultimately, we want to create a map of the study area which shows us different "landcovers" (materials). This can broadly be achieved in three steps:
1. Generate labeled landcover data
2. Train a model using labeled data
3. Validate the model
## 1. Generating Labeled Data
A vital step in any machine learning workflow is the generation of labeled data, which we will use to train a model to differentiate between different types of land cover and later to test the model's accuracy. By looking around the study area, we can get a sense of the different land cover classes that we might encounter:
1. Agricultural Land
2. Urban Areas
3. Oil Contamination
Naturally we could subdivide each of these into sub-categories, and there are probably other classes we haven't included that may be present in the study area. The choice of classes is partly informed by the nature of the task at hand. In theory, the most efficient number of classes for this task would be two: oil, and everything else. The problem is that the "everything else" category would be pretty noisy since it would include a wide range of materials, making it harder to distinguish this from oil. In practice, a visual inspection of major landcover classes in the study area is a quick-and-dirty way of getting at roughly the right number of classes. This is also an iterative process: you can start with a set of labeled data, look at the model results, and adjust your sampling accordingly. More on this later.
The main landcover class we're interested in is, of course, oil. Some oil contamination is readily visible from the high resolution satellite basemap; rivers of oil flow from the leaking [Ger Zero refinery](https://zoom.earth/#view=36.947921,42.02871,16z/overlays=heat,labels:off,crosshair). We can draw polygons around the oil contamination like so:
![](images/ger_zero.png)
The same process is applied to agricultural land and urban areas. In general, you want to make sure that you're sampling from all across the study area. I've generated between four to 10 polygons per landcover class in different places. We're now left with a featureCollection composed of polygons for each class. I've named them `oil`, `agriculture`, and `urban`.
However, I don't just want to use all of the pixels contained in these polygons for training. There are several reasons for this. First, it would likely lead to [overfitting](https://en.wikipedia.org/wiki/Overfitting). Second, there are probably over a million pixels between all of the polygons, which would slow things down unnecessarily. Third, I haven't drawn the polygons to be equal sizes across classes, so I could end up with way more points from one class compared to another. It's OK to have some imbalance between classes, but you don't want it to be extreme.
As such, the next step involves taking random samples of points from *within* these polygons. I do so using the `randomPoints` function:
```js
var oil_points=ee.FeatureCollection.randomPoints(oil, 3000).map(function(i){
return i.set({'class': 0})})
var urban_points=ee.FeatureCollection.randomPoints(urban, 1000).map(function(i){
return i.set({'class': 1})})
var agriculture_points=ee.FeatureCollection.randomPoints(agriculture, 2000).map(function(i){
return i.set({'class': 2})})
```
In the first line, I create a new featureCollection called `oil_points` which contains 3000 points sampled from the polygons in the `oil` featureCollection. I then map through each of these points, and set a property called "class" equal to 0. I do the same for the urban and agricultural areas, setting the "class" property of these featureCollections to 1 and 2, respectively. Ultimately, our model will output a raster in which each pixel will contain one of these three values. A value of 0 in the output will represent the model predicting that that pixel is oil, based on the training data; a value of 1 would indicate predicted urban land cover, and 2 predicted agricultural landcover.
Now we want to create one feature collection called "sample", which will contain all three sets of points.
```js
var sample=ee.FeatureCollection([oil_points,
urban_points,
agriculture_points
])
.flatten()
.randomColumn();
```
We've also assigned a property called "random" using the `randomColumn` function. This lets us split our featureCollection into two: one used for training the model, and one used for validation. We'll use a 70-30 split.
```js
var split=0.7
var training_sample = sample.filter(ee.Filter.lt('random', split));
var validation_sample = sample.filter(ee.Filter.gte('random', split));
```
## 2. Training a Model
Having generated labeled training and testing data, we now want to teach an algorithm to associate the pixels in those areas (in particular, their spectral profiles) with a specific landcover class.
The list of points we generated in the previous step contain a label (0: oil, 1: urban, 2: agriculture). However, they do not yet contain any information about the spectral profile of the Sentinel-2 image. The `sampleRegions` function lets us assign the band values from an image as properties to our feature collection. We do this for both the training sample and the validation sample.
```js
var training = image.sampleRegions({
collection: training_sample,
properties: ['class'],
scale: 10,
});
var validation = image.sampleRegions({
collection: validation_sample,
properties: ['class'],
scale: 10
});
```
Each point in the featureCollections above will contain a property denoting each Sentinel-2 band's value at that location, as well as the property denoting the class label.
Now we're ready to train the model. We'll be using a [Random Forest](https://en.wikipedia.org/wiki/Random_forest) classifier, which basically works by trying to separate your data into the specified classes by setting lots of thresholds in your input properties (in our case, Sentinel-2 band values). It's a versatile and widely-used model.
We first call a random forest classifier with 500 trees. More trees usually yields higher accuracy, though there are diminishing returns. Too many trees will result in your computation timing out. We then train the model using the `train` function, which we supply with the training data as well as the name of the property that contains our class labels ("class").
```js
var model = ee.Classifier.smileRandomForest(500)
.train(training, 'class');
```
The trained model now associates Sentinel-2 band values with one of three landcover classes. We can now feed the model pixels it has never seen before, and it will use what it now knows about the spectral profiles of the different classes to predict the class of the new pixel.
```js
var prediction = image.classify(model)
```
`prediction` is now a raster which contains one of three values (0: oil, 1: urban, 2: agriculture). We're only interested in oil, so let's isolate the regions in this raster that have a value of 0, and add them in red to the map:
```js
var oil_prediction=prediction.updateMask(prediction.eq(0))
Map.addLayer(oil_prediction, {palette:'red'}, 'Predicted Oil Conamination')
```
![](images/ger_zero_pred.png)
## 3. Validation
The image above should look somewhat familiar. It's Ger Zero, where we trained part of our model. We can see in red the areas which the model predicts to be oil pollution. These largely align with the areas that we can see as being contaminated based on the high resolution basemap. It's not perfect, but it's pretty good.
Let's scroll to another area, far from where the model was trained.
![](images/small_refinery.png)
This image shows two clusters of makeshift refineries which were identified by the model. This is good, though we can only get so far by visually inspecting the output from our model. To get a better sense of our model's performance, we can use the validation data that we generated previously. Remember, these are labeled points which our model was not trained on, and has never seen before.
We'll take the `validation` featureCollection containing our labeled points, and have our model classify it.
```js
var validated = validation.classify(model);
```
Now the `validated` variable is a featureCollection which contains both manual labels and predicted labels from our model. We can compare the manual labels to the predicted output to get a sense of how well our model is performing. This is called a Confusion Matrix (or an Error Matrix):
```js
var testAccuracy = validated.errorMatrix('class', 'classification');
print('Confusion Matrix ', testAccuracy);
```
| | | | *Labels* | |
|:------------:|:---------------:|:-------:|:---------:|:---------------:|
| | | **Oil** | **Urban** | **Agriculture** |
| | **Oil** | 876 | 1 | 5 |
| *Prediction* | **Urban** | 0 | 168 | 8 |
| | **Agriculture** | 1 | 4 | 514 |
Now, we can see that of the 877 points that were labeled "oil", only one was falsely predicted to be agricultural land. The model also falsely predicted as oil one point that was labeled urban, and five points that were labeled agriculture. Not bad. We can get a sense of the model's overall accuracy using the `accuracy` function on the confusion matrix:
```js
print('Validation overall accuracy: ', testAccuracy.accuracy())
```
This tells us that the overall accuracy of our model is around 98%. However, we shouldn't take this estimate at face value. There are a number of complicated reasons ([spatial autocorrelation](https://www.sciencedirect.com/topics/computer-science/spatial-autocorrelation#:~:text=Spatial%20autocorrelation%20is%20the%20term,together%20to%20have%20similar%20values.) in the training data, for example) why this figure is probably inflated. If we were submitting this analysis to a peer-reviewed journal, we'd take great care in addressing this, but for our purposes we can use the accuracy statistics to guide our analysis and get a rough sense of how well the model is performing.
This model isn't perfect; it often mis-classifies the shorelines of lakes as oil, or certain parts of urban areas. As previously mentioned, training a model is often an iterative process. At this stage, if your accuracy is not as high as you'd like it to be, you can use the output to figure out how to tweak the model. For example, you may observe that your model is confusing urban areas with oil spills. You can draw a polygon over the erroneous area, label it urban landcover and retrain the model thereby hopefully improving accuracy. We could further refine our model in this way.
# Communicating the Results
Now that we've got a model that can identify oil from multispectral satellite imagery fairly well, we can set about making our results accessible.
One of the things we're particularly interested in is the distribution of small refineries. The way we're currently visualizing the prediction (the raster output from the model where predicted oil is shown in red and everything else is transparent) makes it hard to see these small refineries when we zoom out:
![](images/big_red.png)
We can convert our raster into a series of points using the `reduceToVectors` function. In essence, this takes homogenous regions of an image (e.g., an area predicted to be oil surrounded by an area not predicted to be oil) and converts it into a point:
```js
var vectors = oil_prediction.reduceToVectors({
geometry: AOI,
scale: 10,
geometryType: 'centroid',
eightConnected: true,
labelProperty: 'classification',
maxPixels:1653602926
}).filterBounds(AOI)
Map.addLayer(vectors.style({color: 'black', fillColor: '#00f2ff', pointSize:5}),{},'Oil Contamination Points',false)
```
Now the distribution of small refineries is much more easily visible as blue dots:
![](images/points.png)
If we zoom out even further, we can see clusters of points that correspond to areas of high oil production. Using geolocated photographs, we can roughly ground-truth the model output:
![](images/UNEP.PNG)

578
chapters/C3_Blast.qmd Normal file
View File

@@ -0,0 +1,578 @@
---
title: Blast Damage Assessment
---
On August 4th, 2020, a warehouse full of fertilizer in the port of Beirut exploded. The blast killed over 200 people, injured thousands, and caused widespread damage to the city:
<div style="padding-bottom:56.25%; position:relative; display:block; width: 100%">
<iframe width="100%" height="100%"
src="https://www.youtube.com/embed/LNDhIGR-83w?start=29"
frameborder="0" allowfullscreen="" style="position:absolute; top:0; left: 0">
</iframe>
</div>
Assessing blast damage is a common task for open source investigators, and satellite imagery analysis is one of the best tools we have at our disposal to analyze this sort of phenomenon. [NASA](https://earthobservatory.nasa.gov/images/147098/scientists-map-beirut-blast-damage) used Sentinel-1 imagery in the aftermath of the explosion to generate an estimated damage map. They explain that Sentinel-1 Synthetic Aperture Radar (SAR) imagery is good for this sort of task:
"SAR instruments send pulses of microwaves toward Earth’s surface and listen for the reflections of those waves. The radar waves can penetrate cloud cover, vegetation, and the dark of night to detect changes that might not show up in visible light imagery. When Earth’s crust moves due to an earthquake, when dry land is suddenly covered by flood water, or when buildings have been damaged or toppled, the amplitude and phase of radar wave reflections changes in those areas and indicates to the satellite that something on the ground has changed."
The NASA team produced this estimate the day after the explosion, which is very impressive. However, due to the quick turnaround, they were pretty light on the description of their methodology (they didn’t provide any code, or even explicitly say how exactly they generated the estimate), making their analysis hard to replicate. They also failed to validate their results, which is a critical step in any analysis.
In this case study we'll be developing our own change detection algorithm from scratch, applying to Sentinel-1 imagery of Beirut before and after the blast, and validating our results using building footprints and U.N. damage estimates as the ground truth. Below is the final result of the analysis, which shows building footprints colored according to the predicted level of damage they sustained from the blast:
:::{.column-page}
<iframe src='https://ollielballinger.users.earthengine.app/view/beirutsar' width='100%' height='700px'></iframe>
:::
## Change Detection
There are many ways to detect change between two images. The simplest way would be to take an image taken before an event, and subtract it from an image taken after. This is a good way to get a general sense of where change has occurred, but if you only use two images (one before an event and another after), it would be difficult to differentiate between areas that have changed as a result of the event in question, and areas that have changed for other reasons. Things in Beirut (and cities in general) are constantly changing: construction, cars/planes/ships moving, vegetation growing, etc. So we wouldn't know if the change we're seeing is a result of the explosion or whether that area is generally prone to change. We can overcome this by comparing a bunch of pre-event images to a bunch of post-event images. This way we can see if the change we're seeing is consistent across all of the images. If it is, then we can be fairly confident that the change is a result of the event in question. The mean is simply the sum of all the values ($x_i$) in a set divided by the number of values ($n$):
$$\large \overline{x} = \frac{1}{n} \sum_{i=1}^n x_i$$
But if we just take the average pixel value before and subtract the average pixel value after, we're not accounting for the *variability* of that pixel's values. For example, if we have a pixel that has had an average value of 1 for the month before the event, and a value of 2 in the month after the event, the difference is 1. If that pixel's value is extremely consistent (it never varies by more than 0.1), such a change would be very significant. But if that pixel's value is very variable (it varies by 2 or even 3 on a regular basis), then the change is not significant. So we need to account for the variability of the pixel's values using the standard deviation. It is calculated as the square root of the variance, which is the average of the squared differences from the mean:
$$\large s = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (x_i - \overline{x})^2}$$
With just the mean and the standard deviations of two sets of numbers, we can use a statistical test to determine whether the change in means is significant. The simplest way to do this is to use a pixelwise t-test, which is basically just a signal-to-noise ratio: it calculates the difference between two sample means (signal), and divides it by the standard deviations of both samples (noise). In this case, the two samples are the pre- and post-event images. The t-test is applied to each pixel in the image, allowing us to determine whether the change is statistically significant. Given two groups, $x_1$ before the event, and $x_2$ after the event, the $t$ statistic is calculated as:
$$ \Large t = {\frac{\overline{x_1}-\overline{x_2}} {\sqrt{\frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}}}} $$
Where:
* $\overline{x}$: Sample Mean
* $s^2$: Sample Standard Deviation
* $n$: Number of observations
This procedure gives us a number called a t-value, which is a measure of how many standard deviations the difference between the two means is. We're not going to get into the details here, but a rule of thumb is that if the t-value is greater than 2, then the difference between the two means is significant. If the t-value is less than 2, then the difference is not significant. We're going to calculate the t-value for each pixel in the image to determine whether that pixel has changed significantly following the event in question. You don't need to know the details of the t-test to understand the results (but hopefully you've got an intuition for what it's doing). If you're interested in learning more about statistical tests of this sort, I teach a course on Data Science at the University College London, and have made all of the lectures and courseware open-source. The T-test lecture is here.
## Implementing a t-test in Earth Engine
Now lets go about implementing this in Earth Engine. We'll start by centering the map on the port of Beirut, and setting the map to satellite view, and defining an area of interest (AOI) as a 3km buffer around the port:
```js
Map.setCenter(35.51898, 33.90153, 15);
Map.setOptions("satellite");
var aoi = ee.Geometry.Point(35.51898, 33.90153).buffer(3000);
```
Next, let's define a function in earth engine that will perform the T-Test. The block of code below defines a function to implement a t-test for every pixel in a set of images. The function will be called 'ttest', and takes four arguments:
* s1: the image collection
* shock: the date of the event
* pre_interval: the number of months before the event
* post_interval: the number of months after the event
The function will return a t-value image, which we can use to determine whether a pixel has changed significantly following the event in question.
```js
function ttest(s1, shock, pre_interval, post_interval) {
// Convert the shock date to a date object
var shock = ee.Date(shock);
// Filter the image collection to the pre-event period
var pre = s1.filterDate(
shock.advance(ee.Number(pre_interval).multiply(-1), "month"),
shock
);
// Filter the image collection to the post-event period
var post = s1.filterDate(shock, shock.advance(post_interval, "month"));
// Calculate the mean, standard deviation, and number of images for the pre-event period
var pre_mean = pre.mean();
var pre_sd = pre.reduce(ee.Reducer.stdDev());
var pre_n = ee.Number(pre.filterBounds(aoi).size());
// Calculate the mean, standard deviation, and number of images for the pre-event period
var post_mean = post.mean();
var post_sd = post.reduce(ee.Reducer.stdDev());
var post_n = ee.Number(post.filterBounds(aoi).size());
// Calculate the pooled standard deviation
var pooled_sd = pre_sd
.multiply(pre_sd)
.multiply(pre_n.subtract(1))
.add(post_sd.multiply(post_sd).multiply(post_n.subtract(1)))
.divide(pre_n.add(post_n).subtract(2))
.sqrt();
// Calculate the denominator of the t-test
var denom = pooled_sd.multiply(
ee.Number(1).divide(pre_n).add(ee.Number(1).divide(post_n)).sqrt()
);
// Calculate the Degrees of Freedom, which is the number of observations minus 2
var df = pre_n.add(post_n).subtract(2);
print("Number of Images: ", df);
// Calculate the t-test using the:
// mean of the pre-event period,
// the mean of the post-event period,
// and the pooled standard deviation
var change = post_mean
.abs()
.subtract(pre_mean.abs())
.divide(denom)
.abs()
.subtract(2);
// return the t-values for each pixel
return change;
}
```
An important detail in the code above is that we've actually tweaked the t-test slightly, in two ways.
First, the algorithm above returns tha *absolute* value of t (i.e. the absolute value of the difference between the two means). This is because we're interested in whether the pixel has changed at all, not whether it's changed in a particular direction. Second, we've subtracted 2 from the t-value.
The t-value is a measure of how many standard deviations the difference between the two means is. Generally speaking, if the t-value is greater than 2, then the difference between the two means is considered statistically significant. 2 is a fairly arbitrary cutoff, but it's the most commonly used one since it corresponds to the 95% confidence interval (i.e., there’s less than a 5% chance of observing a difference that big due to random chance). Now we've got a function that can take an image collection and return a t-value image, where a value greater than 0 corresponds to a statistically significant change between the pre-event and post-event periods.
## Filtering the Sentinel-1 Imagery
We can't just blindly apply this algorithm to the entire image collection, because the image collection contains images from both ascending and descending orbits. We need to filter the image collection to the ascending and descending orbits, and then calculate the t-value for each orbit separately: this is because the satellite is viewing the scene from a completely different angle when it's in ascending and descending orbits, which will generate a lot of noise in our data. In fact, even when the satellite is either ascending or descending, we can have multiple images of the same place taken from slightly different orbital tracks because these overlap (see [this visualization of orbits](/ch1#orbits)). We need to filter the image collection to the relative orbit number that is most common within the image collection. For that, we define a new function called 'filter_s1', which takes a single argument: the path (either 'ASCENDING' or 'DESCENDING').
```js
function filter_s1(path) {
// Filter the image collection to the ascending or descending orbit
var s1 = ee
.ImageCollection("COPERNICUS/S1_GRD")
.filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VH"))
.filter(ee.Filter.eq("instrumentMode", "IW"))
.filter(ee.Filter.eq("orbitProperties_pass", path))
.filterBounds(aoi)
.select("VH");
// Find the most common relative orbit number
var orbit = s1
.aggregate_array("relativeOrbitNumber_start")
.reduce(ee.Reducer.mode());
// Filter the image collection to the most common relative orbit number
var s1 = s1.filter(ee.Filter.eq("relativeOrbitNumber_start", orbit));
// Calculate the t-test for the filtered image collection using the function we defined earlier
var change = ttest(s1, "2020-08-04", 12, 2);
// Return the t-values for each pixel
return change;
}
```
You'll notice that we've called the ttest function we defined earlier with four arguments:
* s1: the Sentinel-1 image collection filtered to the ascending or descending orbit
* 2020-08-04: the date of the explosion
* 12: the number of months before the explosion to use for the pre-event period; I'm choosing to include the full year prior to the explosion to get a good baseline
* 2: the number of months after the explosion to use for the post-event period; I'm including 2 months after the explosion. Much less than that, and we risk missing the effects of the explosion. Much more than that, and we risk including the effects other changes that happened after the explosion, including the reconstruction effort.
Now we want to apply this function to the image collection twice (once for each orbit) and then combine the two images into a single image. After that, we can clip it to the area of interest and display it on the map:
```js
// Call the filter_s1 function twice, once for each orbit, and then combine the two images into a single image
var composite = ee
.ImageCollection([filter_s1("ASCENDING"), filter_s1("DESCENDING")])
.mean()
.clip(aoi);
// Define a color palette
var palette = ["440154", "3b528b", "21918c", "5ec962", "fde725"];
// Add the composite to the map
Map.addLayer(
composite,
{ min: 0, max: 4, opacity: 0.8, palette: palette },
"change"
);
```
The visualization parameters correspond to the statistical significance of the change in pixel values. Using the Viridis color palette which ranges from purple to yellow, dark purple pixels indicate no significant change, and yellow pixels indicate a significant change with with 95% confidence. The brighter the yellow of a pixel, the more significant the change.
![Pixelwise T-Test, 2020](../images/beirut/beirut_change_2020.jpg)
This seems to be working quite well; but remember, ports are generally prone to change. The t-test is accounting for this by calculating each pixel's variance over the entire time period, but it's still possible that the change we're seeing is due to the port rather than the explosion. To test this, we can run the same algorithm on the same area, using the same date cutoff (August 4th), but in a different year; I've chosen 2018. This is what's known as a placebo test: if it's still showing loads of statistically significant change around the cutoff, our algorithm is probably picking up on port activity rather than the explosion.
![Pixelwise T-Test, 2018](../images/beirut/beirut_change_2018.jpg)
Compared to the 2020 image, there's a lot less yellow (significant change). That being said, there are a few yellow areas. This could be due to a number of reasons: ships coming and going, cranes moving, and containers being loaded and unloaded would all register in the change detection algorithm. There are also a number of yellow specks throughout the city, which is also to be expected since cities are also generally in a state of flux. Construction, demolition, and even the growth of vegetation can all be detected by the algorithm.
However, the scale and quantity of the change is nowhere near what it was for the 2020 image. This is a good sign that the algorithm is detecting change resulting from the explosion.
## Validation
Great. We've developed our very own change detection algorithm in earth engine, applied it to the Beirut explosion, and it seems to be working after checking with a basic placebo test. But how do we know that it's correctly predicting the *extent* of the damage, and not wildly over/underestimating?
Given that this was a few years ago, we have the benefit of hindsight. In particular, the United Nations and the Municipality of Beirut have [published a report](https://unhabitat.org/sites/default/files/2020/10/municipality_of_beirut_-_beirut_explosion_rapid_assessment_report.pdf) on the damage caused by the explosion. This report includes estimates of the number of buildings damaged or destroyed, as well as the number of people displaced. The report states that approximately 10,000 buildings were damaged within a 3 kilometre radius of the port. If our algorithm suggests that only 1,000 buildings were damaged, it's undershooting. If it suggests that 100,000 buildings were damaged, it's overshooting.
Using building footprint data and the t-test image we just generated, we can createe an estimate of the number of damaged buildings according to our model. First, we want to generate a thresholded image, where pixels with a value greater than 0 are set to 1, and all other pixels are set to 0. We can then use this mask to reduce the building footprints to a single value for each building, where the value is the mean of the t-test image within the footprint. If the mean value is greater than 0, the building is damaged. If it's less than 0, the building is not damaged.
```js
// Create a mask of the t-test image, where pixels with a value greater than 0 are set to 1, and all other pixels are set to 0
var threshold = composite.updateMask(composite.gt(0));
// Load the building footprints
var buildings = ee
.FeatureCollection("projects/sat-io/open-datasets/MSBuildings/Lebanon")
.filterBounds(aoi);
// Calculate the mean value of the t-test image within each building footprint
var damaged_buildings = threshold.reduceRegions({
collection: buildings,
reducer: ee.Reducer.mean(),
scale: 1,
});
// Print the number of buildings with a mean value greater than 0
// i.e., those displaying statistically significant change
print(damaged_buildings.filter(ee.Filter.gt("mean", 0)).size());
```
The result is 9,256, which is pretty damn close to 10,000. We can also visualize the building footprints on the map, colored according the mean value of the t-test image within the footprint, where:
* Blue = no damage
* Green = low damage
* Yellow/Orange = medium damage
* Red = high levels of damage
```js
// Create an empty image
var empty = ee.Image().byte();
// Paint the building footprints onto the empty image
var outline = empty.paint({
featureCollection: damaged_buildings,
color: "mean",
width: 5,
});
// Define a color palette
var building_palette = [
"0034f5",
"1e7d83",
"4da910",
"b3c120",
"fcc228",
"ff8410",
"fd3000",
];
// Add the image to the map
Map.addLayer(
outline,
{ palette: building_palette, min: 0, max: 2 },
"Damaged Buildings"
);
```
The result naturally resembles the underlying t-test image, with high levels of damage concentrated around the port, and progressively decreasing damage with distance:
![Building Footprints colored according to estimated blast damage](../images/beirut/beirut_footprints.jpg)
To get a better sense of how these predicitions correspond to actual damage, we can zoom in and turn on the Google satellite basemap, which has imagery taken just after the explosion; you can still see capsized boats in the port. Zooming in to the epicenter, we can see several warehouses that were effectively vaporized. Our change detection algorithm picks up on a high degree of change, as indicated by the red outlines of the building footprints:
![Predicted damage and optical satellite imagery in the Port of Beirut, August 2020](../images/beirut/beirut_footprints_port.jpg)
This is pretty low-hanging fruit. Let's look at a different area, around 1.3km east from the epicenter with a mix of warehouses and residential buildings:
![Area east of the port: 35.533194, 33.9024](../images/beirut/beirut_footprints_zoomed.jpg)
Here, there's greater variation in the predictions. I've highlighted three areas.
In Area A, we see a warehouse with a highly deformed roof; panels of corrugated iron are missing, and much of the roof is warped. The building footprint for this warehouse is red, suggesting that our algorithm correctly predicts a significant amount of blast damage.
In Area B, we see a medium-rise building. If you look closely at the southern edge of the building, you'll see the siding has been completely torn off and is laying on the sidewalk. The building footprint is orange, suggesting a medium amount of change. We may be underestimating a bit here.
In Area C, there are a bunch of high rise buildings clustered together. The building footprints are all blue, suggesting little to no damage. This is a bit of a surprise given how damaged areas A and B are. If you squint at the satellite image, it is indeed hard to tell if these buildings are damaged because we're looking at them from the top down, when much of the damage (e.g., the windows being blown out) would only be visible from the side. Indeed, our own estimate of the number of damaged buildings based on the algorithm we developed is about 8% shy of the U.N.'s estimate. This may be why.
## Conclusion
In this practical example, we created a custom change detection algorithm that conducts a pixelwise t-test to detect change resulting from the 2020 explosion in the port of Beirut. By defining our own functions to do most of this analysis, we can apply the same workflow quite easily to a different context by simply moving the AOI and inputting the date of the shock. A placebo test showed that it's not just detecting general change in the area, but specifically change resulting from the explosion: when we keep everything the same but change the year of the shock, we see very little significant change being detected. Finally, by joining the predicted damage map to building footprints, we come up with an estimate of 9,256 damaged buildings, which is pretty close to the U.N.'s estimate of 10,000. That concludes the portion of this case study that deals with Earth Engine, but if you're interested in learning more about why we're coming up a bit short on the damage estimate (and some different ways of looking at the problem), read on.
## Extension: Satellite Imagery and its Limits
Though satellite imagery analysis is undoubtedly one of the best tools we have at our disposal to analyze this sort of phenomenon, it appears to systematically underestimate the extent of damage in Beirut. I outline an alternative approach using Open Street Map data to create a 3D model of Beirut and the explosion to analyze directional blast damage. Again, we're now leaving Earth Engine and moving to Blender, so if you're not interested in that feel free to skip ahead to the next case study.
Below is one of the most viewed videos of the explosion:
<blockquote class="twitter-tweet"><p lang="en" dir="ltr">Stunning video shows explosions just minutes ago at Beirut port <a href="https://t.co/ZjltF0VcTr">pic.twitter.com/ZjltF0VcTr</a></p>&mdash; Borzou Daragahi 🖊🗒 (@borzou) <a href="https://twitter.com/borzou/status/1290675854767513600?ref_src=twsrc%5Etfw">August 4, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>
Geolocating this video was pretty simple thanks to the Greek Orthodox church (highlighted in green below) and the road leading to it (highlighted in blue). The red box indicates the likely location (33.889061, 35.515909) from which the person was filming:
![](../images/beirut/IMG_2.png)
The video shows heavy damage being sustained by areas well outside the zones classified as damaged in the maps above (both my own and NASA's). Indeed, substantial damage was reported several kilometers away.
Why are satellite images underestimating damage in Beirut? Satellite images are taken from above, and are two-dimensional. Much of the damage caused by the blast, however, was directional; the pressure wave hit the sides of buildings, as shown in this diagram from a FEMA manual:
![](../images/beirut/IMG_3.png)
Areas close to the explosion suffered so much damage that it could be seen from above, but even if an apartment building had all of its windows blown out, this would not necessarily be visible in a top-down view. Even for radar, which does technically collect data in three dimensions, the angle problem remains; a high resolution radar might be able to tell you how tall an apartment complex is, but it won't give you a clear image of all sides. Case in point: the NASA damage map was created using Sentinel-1 SAR data. In a nutshell, damage assessment in this case is a three-dimensional problem, and remote sensing is a two-dimensional solution.
### Creating a 3D model of Beirut
To create a more accurate rendering of directional blast damage, three dimensional data are required. Data from Open Street Maps (OSM) contains information on both the "footprints" (i.e., the location and shape) as well as the height of buildings, which is enough to create a three dimensional model of Beirut. 3D rendering was done in Blender using the Blender-OSM add-on to import a satellite basemap, terrain raster, and OSM data.
Geolocated videos of the blast can be used to verify and adjust the model. Below is a side-by-side comparison of the twitter video and a 3D rendition of OSM data:
![](../images/beirut/IMG_4.png)
Some slight adjustments to the raw OSM data were made to achieve the image on the right. The building footprints are generally very accurate and comprehensive in coverage, but the building height data does occasionally have to be adjusted manually. A simple and reliable way of doing this is to look at the shadows cast by the building on the satellite base map and scale accordingly. I also added a rough texture to the buildings to help differentiate them, and added the domed roof of the Greek Orthodox church for reference.
For good measure, a second video is geolocated following the same procedure:
<blockquote class="twitter-tweet"><p lang="en" dir="ltr">Another view of the explosions in Beirut <a href="https://t.co/efT5VlpMkj">pic.twitter.com/efT5VlpMkj</a></p>&mdash; Borzou Daragahi 🖊🗒 (@borzou) <a href="https://twitter.com/borzou/status/1290678580897251330?ref_src=twsrc%5Etfw">August 4, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>
The second pier (highlighted in green) and the angle (in blue) serve as references:
![](../images/beirut/IMG_5.png)
The video was taken from the rooftop of a Japanese restaurant called Clap Beirut (in red above). This is confirmed by a picture of the rooftop bar on google images, which matches the bar that can be seen at 0:02 in the twitter video. Below is a comparison of the video view and the 3D OSM model:
![](../images/beirut/IMG_6.png)
Though somewhat grainy, the basemap on the OSM rendering shows the same parking lot in the foreground, the second pier, and the same two buildings highlighted in yellow. Having created a 3D model of Beirut using OSM data, we can now simulate how the explosion would interact with the cityscape.
### Using a Viewshed Analysis to Assess Blast Exposure
As the pressure wave moved through the Beirut, some buildings bore the full force of the explosion, while others were partially shielded by taller structures. A viewshed analysis can be conducted to identify surfaces that were directly exposed to the explosion by creating a lighting object at ground zero; areas that are lit up experienced unobstructed exposure to the blast:
![](../images/beirut/GIF_1.gif)
Pressure waves, like sound, are capable of diffraction (bending around small obstructions). To roughly simulate this, the lighting object is gradually raised, allowing the light to pass "around" obstructions. Warehouses on the Eastern side of the docks, as well as the first row of apartment buildings facing the docks are immediately affected. As the lighting object rises above the warehouse, more areas suffer direct exposure.
Using two lighting objects-- a red one at 10 meters and a blue one at 20 meters above the warehouse at ground zero-- the intensity of the blast in different areas is highlighted; red areas suffered direct exposure, blue areas suffered partially obstructed exposure, and black areas were indirectly exposed.
![](../images/beirut/IMG_7.png)
In the immediate vicinity of the explosion the large "L" shaped building (Lebanon's strategic grain reserve) is bright red, and was barely left standing. It absorbed a large amount of the blast, shielding areas behind it and thereby casting a long blue shadow to the West. If one refers back to the satellite damage maps above, there appears to be significantly less damage in the area just West of ("behind") the grain silo, roughly corresponding to the blue shadow above. While these areas were still heavily damaged, they seem to have suffered less damage than areas of equal distance to the East.
### Accounting for Diffraction
The viewshed analysis tells us which sides of a building are exposed to the blast, but it's a pretty rough approximation of the way the pressure wave would respond to obstacles in its path. As previously mentioned, pressure waves behave much like sound waves or waves in water: they bounce off of objects, move around obstructions, and gradually fade.
To get a more precise idea of the way in which the blast interacted with the urban environment, we can model the blast as an actual wave using the "dynamic wave" feature in Blender. This effectively involves creating a two-dimensional plane, telling it to behave like water, and simulating an object being dropped into the water. By putting an obstruction in this plane, we can see how the wave responds to it. As an example, the grain silo has been isolated below:
![](../images/beirut/GIF_2.gif)
As the blast hits the side of the silo, it is reflected. Two large waves can be seen traveling to the right: the initial blast wave, and the reflection from the silo which rivals the initial wave in magnitude. To the left, the wave travels around the silo but is significantly weakened.
Broadening the focus and adding the rest of the OSM data back in, we can observe how the pressure wave interacted with buildings on the waterfront:
![](../images/beirut/GIF_3.gif)
The warehouses on the docks were omitted to emphasize the interaction between the pressure wave and the waterfront buildings; their light metal structure and low height means they would have caused little reflection anyway. The general pattern of the dynamic wave is consistent with the viewshed, but adds a layer of detail. The blast is reflected off of the silo towards the East, leading to a double hit. Though the wave still moves around the silo to the West, the pressure is diminished. Once the wave hits the highrises, the pattern becomes noisy as the wave both presses forward into the mainland and is reflected back towards the pier.
### Modeling the Pressure Wave
Now that we've accounted for the directionality of the blast and the influence of buildings, we can model the pressure wave itself. An expanding sphere centered at ground zero is used to model the progression of the pressure wave through the city. To get a visual sense of the blast's force, the color of the sphere will be a function of the pressure exerted by the pressure wave.
The pressure exerted by the explosion in kilopascals (kPa) at various distances can be calculated using the U.S. Department of Defemse's Blast Effects Computer, which allows users to input variables such as the TNT equivalent of the ordnance, storage method, and elevation. Though there are several estimates, the blast was likely equivalent to around 300 tons of TNT. The direct "incident pressure" of the pressure wave is shown in blue. However, pressure waves from explosions that occur on the ground are reflected upwards, amplifying the total pressure exerted by the blast. This "reflected pressure" is shown in orange:
<iframe title="Blast Overpressure and Distance " aria-label="Interactive line chart" id="datawrapper-chart-J1Pb1" src="https://datawrapper.dwcdn.net/J1Pb1/1/" scrolling="no" frameborder="0" style="width: 0; min-width: 100% !important; border: none;" height="400"></iframe><script type="text/javascript">!function(){"use strict";window.addEventListener("message",(function(a){if(void 0!==a.data["datawrapper-height"])for(var e in a.data["datawrapper-height"]){var t=document.getElementById("datawrapper-chart-"+e)||document.querySelector("iframe[src*='"+e+"']");t&&(t.style.height=a.data["datawrapper-height"][e]+"px")}}))}();
</script>
For reference, 137 kPa results in 99% fatalities, 68 kPa is enough to cause structural damage to most buildings, and 20 kPa results in serious injuries. 1-6 kPa is enough to break an average window. At 1km, the reflected pressure of the blast (18 kPa) was still enough to seriously injure. Precisely calculating the force exerted by an explosion is exceptionally complicated, however, so these numbers should be treated as rough estimates. Further analysis of the damage caused by the blast can be derived from the UN's Explosion Consequences Analysis calculator which provides distance values for different types of damage and injuries.
Linking the values in this graph to the color of the pressure wave sphere provides a visual representation of the blast's force as it expands. An RGB color scale corresponds to the blast's overpressure at three threshold values.
![](../images/beirut/beirut.gif)
By keeping the lighting object from the viewshed analysis and placing it within the expanding sphere of the pressure wave, we combine two key pieces of information: the pressure exerted by the blast (the color of the sphere), and the level of directional exposure (brightness).
Now, referring back to the two geolocated twitter videos from earlier, we can recreate the blast in our 3D model and get some new insights. Below is a side-by-side comparison of the first video and the 3D model:
![](../images/beirut/GIF_5.gif)
Judging by the twitter video alone, it would be very hard to tell the fate of the person filming or the damage caused to the building that they were in. However, the 3D model shows that despite having an unobstructed view of the explosion, the incident pressure of the pressure wave had decreased significantly by the time it reached the viewing point. The blue-green color corresponds to roughly 15 kPa-- enough to injure and break windows, but not enough to cause structural damage to the building.
The second twitter video was taken slightly closer to ground zero, but the view was partially obstructed by the grain silo:
![](../images/beirut/GIF_6.gif)
Though the pressure wave probably exerted more pressure compared to the first angle, the partial obstruction of the grain silo likely tempered the force of the blast.
### Assessing Damage to the Skyline Tower
As a concrete example of how this approach can be used to assess damage (or predict it, if one had the foresight), let us consider the Skyline Tower, pictured below following the explosion:
![](../images/beirut/IMG_8.png)
This partial side view shows two faces of the building, labeled "A" and "B" above. Side A was nearly perpendicular to the blast, and just 600m from ground zero. Based on the previous modeling, the pressure wave exerted roughly 40 kPa on this side of the building. The corner where sides A and B meet, highlighted in green, shows total destruction of windows, removal of most siding panels, and structural damage. The back corner, highlighted in red, shows many windows still intact, indicating that the maximum overpressure on this side of the building likely didn't exceed 10 kPa. In other words, standing on the front balcony would likely have led to serious injury but standing on the back balcony would have been relatively safe.
The animation below shows the Skyline Tower as it is hit by the pressure wave, with sides A and B labeled:
![](../images/beirut/GIF_7.gif)
The bright green color of the pressure wave indicates a strong likelihood of structural damage. Side A can be seen taking a direct hit, while side B is angled slightly away. Despite not being directly exposed to the blast, it likely still took reflective damage from some of the neighboring buildings. Both the incident overpressure indicated by the color of the sphere, as well as the relative brightness of sides A and B both correspond closely to the observed damage taken by the Skyline Tower.
### Further Research
Though satellite imagery analysis is an indispensable tool in disaster response, it has limitations. Urban blast damage in particular is difficult to assess accurately because it is highly directional and much of it cannot be seen from a bird's eye view. Using free and open source tools, an interactive 3D model of an urban explosion can be generated, allowing for a highly detailed investigation of directional blast damage. This can be achieved in three steps:
First, creating a 3D model of the urban area using Blender and Open Street Maps data. Second, conducting a viewshed analysis using lighting objects to gauge levels of unobstructed exposure to the pressure wave. Third, modeling the explosion using geolocated videos of the event and ordnance calculators. For added detail, a dynamic wave analysis can be used to more precisely model how the pressure wave interacts with buildings.
Once properly modeled, the explosion can be viewed from any angle in the city. The viewshed analysis can be calibrated more finely by ground-truthing various damage levels (e.g. broken windows) at different locations. In the absence of an official address registry in Beirut, OSM is already being used by the Lebanese Red Cross (donate here) to conduct neighborhood surveys assessing blast damage. As such, this type of damage analysis can quickly be integrated into relief efforts, adapted to model disasters in different cities, and can even be used to simulate the destructive potential of hypothetical explosions to promote readiness.

465
chapters/C4_Ships.qmd Normal file
View File

@@ -0,0 +1,465 @@
---
title: Ship Detection
---
There's a huge amount of data available on the internet about ship movements, most of which draws on the Automatic Identification System (AIS) which uses radio to broadcast the identity, position, course, speed and other data about ships. [MarineTraffic](https://www.marinetraffic.com/en/ais-api-services), for example, provides an API that allows you to query the location of ships in real time as well as historical vessel tracks and lots of other useful data. Unfortunately most sources of AIS data are paywalled, and AIS can be turned off or manipulated to hide the identity or position of the ship. In fact, most of the stuff we're interested in investigating probably happens when AIS is turned off.
Though ships can hide by turning off their AIS transponders, they can't hide from satellites. In this tutorial, we're going to build an application that uses Synthetic Aperture Radar (SAR) from the European Space Agency's Sentinel-1 satellite to automatically identify ships, regardless of whether they've got their transponders turned on or off. Here's the finished application:
:::{.column-page}
<iframe src='https://ollielballinger.users.earthengine.app/view/shipdetection' width='100%' height='700px'></iframe>
:::
## How it Works
The app has two main panels:
1. A control panel on the left that allows the user to interact with the application
2. A map on the right that displays the results
The control panel has a date slider that allows the user to load imagery from a particular year. Below that is a graph that shows the number of ships detected over time within that year. A slider underneath the graph lets us toggle the sensitivity of the ship detection process. Finally, a button at the bottom lets the user draw their own area of interest on the map, and the app will automatically detect ships within that area.
The map panel visualizes the results of the ship detection process and has three layers. The bottom layer is the Sentinel-1 image that we're using to detect ships; it's blue/purple, and if you zoom in and look cloesly you can see bright specks in the sea, which are ships. When Sentinel-1 sends a pulse of radio waves onto a flat surface like the sea, there is very little to reflect the waves back to the satellite -- they just bounce off into space. A low return signal means we'll see a darker color on our map. But when the radio waves hit a ship they are reflected back to the satellite and generate a higher return signal, and therefore a much brighter color. The second layer on the map displays a bunch of green points; each one of these is a detected ship. The last layer shows the red outline of the area of interest that the user drew on the map. You can zoom in on the map by holding down the command button and scrolling up and down.
When the application is first loaded it is centered on an area just north of the Suez Canal, and is analyzing imagery from 2021. We can see a bunch of green dots in the AOI, which is the main waiting area for ships waiting to transit the canal. It's a bit crowded because it's visualizing all of the ships detected in the entire year. We can display imagery from a single day by clicking on a point in the graph on the left, which you will notice displays a huge spike in the number of ships detected around March.
You might remember that on March 23rd, 2021, the Ever Given -- a 400m long container ship -- got stuck in the Suez Canal. The ship was blocking the canal for six days, and it's estimated that it cost the global economy $400 million per day. If you click on the tip of the spike on March 30th, you can see a backup of around 150 ships waiting for the canal to be cleared. You can also zoom in on a particular date range by scrolling and dragging on the graph. If you zoom in on the spike, you can then select imagery from early April to compare the number of ships in the waiting area after the blockage was cleared. In normal times we can see a regular pattern in the number of ships in the waiting area ranging between 15 and 40 ships.
If you're closely zoomed in to the map and load imagery from different days by clicking on the graph, you can compare the bright spots on the Sentinel image and the green dots. The ship detection process is pretty accurate, and we typically see one green dot per ship. However, you may notice that we occasionally miss a ship. This is because the ship detection process is based on a threshold, and if the ship is too small it may not generate a high enough return signal to be detected. You can increase the sensitivity of the ship detection process by moving the slider below the graph. This will increase the number of ships detected, but it may also increase the number of false positives.
The next section focuses on building this application. After that, we'll have a look at a few different use cases for this sort of maritime surveillance.
# Building the Application
## Setup
The first step is to configure the map and import the necessary datasets. By default, we want the app to be centered on the Suez Canal. Then, we want to import the Digital Surface Model (DSM) from the ALOS World 3D-30 dataset. This dataset provides a 30m resolution elevation model of the Earth which we will use to mask out the land. Finally, we want to import the Sentinel 1 dataset. We will use the VV polarization and the Interferometric Wide (IW) mode. We will also sort the images by date.
```js
// Center the map on the Suez Canal and set map options
Map.setCenter(32.327, 31.4532, 10);
Map.setOptions("Hybrid");
Map.setControlVisibility({ all: false });
// Import the Digital Surface Model (DSM) from the ALOS World 3D-30 dataset
var dem = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2").mean().select("DSM");
// Import the Sentinel 1 dataset
var s1 = ee
.ImageCollection("COPERNICUS/S1_GRD")
.filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VV"))
.filter(ee.Filter.eq("instrumentMode", "IW"))
.sort("system:time_start");
// Define the default area of interest
var suez = ee.Geometry.Polygon([
[
[32.17388584692775, 31.59541178442045],
[32.17388584692775, 31.327159861902278],
[32.4787564523965, 31.327159861902278],
[32.4787564523965, 31.59541178442045],
],
]);
```
Now that we've gotten that out of the way, we can move on to the actual detection of ships.
## Ship Detection
You might expect the automatic identification of ships based on synthetic aperture radar satellite imagery to involve a complex machine learning algorithm or artificial intelligence. In fact, it can be done in one line of code which sets a cutoff. If the return signal is greater than 0, then we have a ship. If it's less than 0, then we don't. Simple as that.
The main analytical function responsible for ship identification is the `getVectors` function shown below. It takes an image as an input and returns a FeatureCollection of points, each corresponding to a ship. The function clips the image to the area of interest, selects the VV polarization, and finally filters out areas where the VV value is smaller than 0. This results in a raster image where the sea is black and the ships are white. We then use the `reduceToVectors` function to convert the raster image to a FeatureCollection of points. The function returns this FeatureCollection, and sets a property called `count` which is the number of ships detected in the image.
```js
function getVectors(img) {
// Get the area of interest from the drawing tools widget.
var aoi = drawingTools.layers().get(0).getEeObject();
// Clip the image to the area of interest
// Select the VV polarization
// Filter areas where the VV value is greater than 0
var cutoff = img.clip(aoi).select("VV").gt(0)
// Convert the raster image to a FeatureCollection of points
var points = cutoff.reduceToVectors({
geometry: aoi,
scale: scaleSlider.getValue(),
geometryType: "centroid",
eightConnected: true,
maxPixels: 1653602926,
});
// Set the number of ships detected in the image as a property called "count"
var count = points.size();
// Set the date of the image as a property called "system:time_start"
var date = ee.Date(img.get("system:time_start"));
return points.set("count", count).set("system:time_start", date);
}
```
The `count` and `system:time_start` properties are used to create the graph of daily ship counts and allow the resulting vector (point) data to interact with the date slider widget. An important detail here is that the "scale" parameter of the `reduceToVectors` function is set to the value of the scale slider widget. This allows the user to adjust the resolution of the ship detection process; a smaller value will allow us to detect smaller ships.
## Visualization
The `viz` function is responsible for displaying the results of the ship detection process. It takes the area of interest, the vector data, and the Sentinel 1 image as inputs. Nothing super complicated here; we're just creating three layers and adding them to the map in order: the underlying Sentinel-1 image raster, the ship vector data in green, and the area of interest outline in red. We're using the `Map.layers().set()` function to replace the existing layers with the new ones, rather than adding new ones each time.
```js
function viz(aoi, vectors, s1Filtered) {
// Create an empty image into which to paint the features, cast to byte.
var empty = ee.Image().byte();
// Paint all the polygon edges with the same number and width, display.
var outline = empty.paint({
featureCollection: aoi,
color: 1,
width: 3,
});
// Create a layer for the area of interest in red
var aoi_layer = ui.Map.Layer(outline, { palette: "red" }, "AOI");
// Create a layer for the vector data in green
var vectorLayer = ui.Map.Layer(
vectors.flatten(),
{ color: "#39ff14" },
"Vectors"
);
// Create a layer for the Sentinel 1 image in false color
var sarLayer = ui.Map.Layer(
s1Filtered,
{ min: [-25, -20, -25], max: [0, 10, 0], opacity: 0.8 },
"SAR"
);
// Add the layers in order
Map.layers().set(0, sarLayer);
Map.layers().set(1, vectorLayer);
Map.layers().set(2, aoi_layer);
}
```
We want a function to handle the visualization because there are two different situations in which we're going to visualize results, and we don’t want to repeat our code. The first situation is when the user draws a new area of interest, moves the date slider or alters the scale. In this case, we want to visualize the results of the ship detection process for the entire year's worth of Sentinel-1 imagery. The second situation is when the user clicks on the chart to analyze a particular day. In this case, we obviously only want to visualize the results of the ship detection process on that day. With this function, we can simply pass the appropriately filtered versions of the Sentinel-1 image and vector data to the function, and it will visualize the results, rather than having to write the same code twice.
## Putting it all together
Having defined a few helper functions to handle the visualization and ship detection process, we can now move on to the main function that will perform the analysis. This will be performed by the `daterangeVectors` function. In a nutshell, it reads the user specified date range from the date slider widget, and filters the Sentinel 1 dataset to only include images within that period. Then, it will loop through each Sentinel-1 image from that year and apply the `getVectors` function to count the number of ships that fall within the area of interest and generate a dataset of points corresponding to detected ships. We'll then use the `viz` function we just defined to visualize all of the ship detections and Sentinel-1 images in the AOI during that year stacked on top of each other. Finally, we'll create a chart based on the number of ships detected per day, and allow the user to click on the chart to visualize the results for a particular day.
```js
var daterangeVectors = function () {
// Get the date range from the date slider widget.
var range = ee.DateRange(
ee.Date(dateSlider.getValue()[0]),
ee.Date(dateSlider.getValue()[1])
);
// Get the area of interest from the drawing tools widget.
var aoi = drawingTools.layers().get(0).getEeObject();
// Hide the user-drawn shape.
drawingTools.layers().get(0).setShown(false);
// Filter the Sentinel 1 dataset to only include images within the date range, and within the area of interest.
var s1Filtered = s1.filterDate(range.start(), range.end()).filterBounds(aoi);
// Count the number of ships in each image using the getVectors function
var vectors = s1Filtered.map(getVectors);
// Use the viz function to visualize the results
viz(aoi, vectors, s1Filtered.max().updateMask(dem.lte(0)));
// Create a chart of the number of ships per day
var chart = ui.Chart.feature
.byFeature({
features: vectors,
xProperty: "system:time_start",
yProperties: ["count"],
})
.setOptions({
title: "Daily Number of Ships in Area of Interest",
vAxis: { title: "Ship Count" },
explorer: { axis: "horizontal" },
lineWidth: 2,
series: "Area of Interest",
});
// Add the chart at a fixed position, so that new charts overwrite older ones.
controlPanel.widgets().set(4, chart);
// Add a click handler to the chart to filter the map by day.
chart.onClick(filterDay);
};
```
There's one function referenced above -- `filterDay`-- that we haven't defined yet. This function is called when the user clicks on the chart to analyze a particular day. It takes the date of the clicked day as an input, filters the Sentinel-1 dataset and vector data accordingly, and uses the `viz` function to display the results for that day.
```js
function filterDay (callback) {
// Get the date of the clicked day
var date = ee.Date(callback);
// Filter the vector data to only include images from that day
var vectorDay = vectors.filterDate(date);
// Filter the Sentinel-1 imagery to only include images from that day
var s1Day = s1.filterDate(date).max().updateMask(dem.lte(0));
// Use the viz function to visualize the results
viz(aoi, vectorDay, s1Day);
};
```
The analytical portion of the application is now complete. Now we have to build a user interface that lets us interact with the application.
## Building a User Interface
There are four main steps in the process of creating the User Interface (UI):
1. Configure the drawing tools that allow the user to draw a polygon on the map.
2. Create some widgets
### Drawing Tools
We eventually want to allow the user to draw a polygon on the map, and count the number of ships that fall within it. In order to do so, we need to set up a few functions related to the drawing tools that allow the user to do this. Among other things, we want to make sure that we're clearing the old geometries so that we're only ever conducting analysis inside the most recent user-drawn polygon, so we'll need to clear the old ones. We also want to specify the type of polygon the user can draw, which for ease will be a rectangle (you could change this to the actual "polygon" type if you wanted to draw more complex geometries).
```js
var drawingTools = Map.drawingTools();
// Remove any existing layers
while (drawingTools.layers().length() > 0) {
var layer = drawingTools.layers().get(0);
drawingTools.layers().remove(layer);
}
// Add a dummy layer to the drawing tools object (the Suez Canal box)
var dummyGeometry = ui.Map.GeometryLayer({
geometries: null,
})
.fromGeometry(suez)
.setShown(false);
// Add the dummy layer to the drawing tools object
drawingTools.layers().add(dummyGeometry);
// Create a function that clears existing geometries and lets the user draw a rectangle
function drawPolygon() {
var layers = drawingTools.layers();
layers.get(0).geometries().remove(layers.get(0).geometries().get(0));
drawingTools.setShape("rectangle");
drawingTools.draw();
}
```
### Widgets
The control panel will eventually contain a few different widgets that allow the user to interact with the application. We'll start by creating a button that allows the user to draw a polygon on the map. We'll also create a slider that allows the user to adjust the size of the ships that are detected (remember, this manipulates the "scale" parameter in the `reduceToVectors` function used in the detection process). The slider will have an accompanying label that tells the user what it does.
```js
// Create a button that allows the user to draw a polygon on the map
var drawButton = ui.Button({
label: "🔺" + " Draw a Polygon",
onClick: drawPolygon,
style: { stretch: "horizontal" },
});
// Create a slider that allows the user to adjust the size of the ships that are detected
var scaleSlider = ui.Slider({
min: 1,
max: 100,
value: 80,
step: 1,
onChange: daterangeVectors,
style: { width: "70%" },
});
// Create a label for the slider
var scaleLabel = ui.Label("Ship Size: ");
// Create a panel that contains the slider and its label
var scalePanel = ui.Panel({
widgets: [scaleLabel, scaleSlider],
style: { stretch: "horizontal" },
layout: ui.Panel.Layout.Flow("horizontal"),
});
```
The last widget we're going to define is the date slider. This widget will trigger the `daterangeVectors` function, which will filter the Sentinel-1 dataset to only include images from the selected year, and then run the detection process on the filtered dataset.
```js
// Specify the start and end dates for the date slider
var start = "2014-01-01";
var now = Date.now();
// Create a date slider that allows the user to select a year
var dateSlider = ui.DateSlider({
value: "2021-03-01",
start: start,
end: now,
period: 365,
onChange: daterangeVectors,
style: { width: "95%" },
});
```
### The Control Panel
Now we're going to assemble all of the widgets we've just defined into one panel, alongside some explanatory text. I'm adding a blank label to the panel as a placeholder for the chart, since it will be re-added to the panel every time the user changes the date on the date slider, the AOI, or the scale.
```js
var controlPanel = ui.Panel({
widgets: [
ui.Label("SAR Ship Detection", {
fontWeight: "bold",
fontSize: "20px",
}),
ui.Label(
"This tool identifies ships using Synthetic Aperture Radar imagery. Use the date slider below to analyze a given year. Click on the graph to show ships on a given day.",
{ whiteSpace: "wrap" }
),
dateSlider,
ui.Label(),
scalePanel,
ui.Label(
"Click the button below and draw a rectangle on the map to count ships in a custom area."
),
drawButton
],
style: {maxWidth: "400px"},
layout: ui.Panel.Layout.flow("vertical", true),
});
```
Once the control panel has been defined, we can add it to
```js
// Add the control panel to the map
ui.root.insert(0,controlPanel);
// Trigger the daterangeVectors function when the user draws a polygon
drawingTools.onDraw(ui.util.debounce(daterangeVectors, 500));
// Run the daterangeVectors function to initialize the map
daterangeVectors();
```
And there we have it. A fully functional, all weather, daytime/nighttime ship detection tool that doesn't rely on AIS data. Let's play around with it.
## Taking it for a spin
### North Korea
In 2020, North Korea implemented one of the most severe COVID-19 lockdowns in the world including a near-total ban on ["all cross-border exchanges, including trade, traffic, and tourism".](https://thediplomat.com/2023/01/north-korea-likely-to-lift-pandemic-border-restrictions-in-2023/). Measures have been so severe that the country appears to have experienced a significant [famine](https://foreignpolicy.com/2022/05/16/kim-north-korea-covid-outbreak-pandemic/). Though there were signs that things have gradually returned to normal, information on North Korea's economy is pretty hard to come by. Ship traffic in and out of the country's largest port, Nampo, is probably a pretty good indicator of the country's economic activity.
But we can't just head on down to Marine Tracker or other services that use AIS data to track ship movements. According to the [U.S. Treasury](https://home.treasury.gov/system/files/126/dprk_vessel_advisory_02232018.pdf), "North Korean-flagged merchant vessels have been known to intentionally disable their AIS transponders to mask their movements. This tactic, whether employed by North Korean-flagged vessels or other vessels involved in trade with North Korea, could conceal the origin or destination of cargo destined for, or originating in, North Korea." They should know -- they're the ones imposing the sanctions that make it illegal to trade with North Korea.
A New York Times [investigation](https://www.nytimes.com/2019/07/16/world/asia/north-korea-luxury-goods-sanctions.html) tracked the maritime voyage of luxury Mercedes cars from Germany to North Korea via the Netherlands, China, Japan, South Korea, and Russia. AIS transponders were turned off at several points throughout this journey, and the investigation had to rely on satellite imagery to fill in the gaps.
Though they used high resolution optical imagery to follow individual ships, we want to identify lots of ships in a large area over a long period. That would get very expensive, and automatic ship detection in optical imagery is relatively difficult. Here's how our SAR tool fares when we draw a box in the bay of Nampo:
![](/images/ships_north_korea.jpg)
Looking at imagery from 2021, we can see ship traffic increasing from nearly zero to around 40 ships per day.
### Ukraine
Odessa is Ukraine's largest port. Following its invasion of Ukraine in February 2022, Russia instituted a naval blockade against Ukrainian ports. The impact of this blockade is clearly visible using the tool we've just built:
![](/images/ships_ukraine.jpg)
The daily number of ships detected in the port of Odessa dropped from 40 to 50 to 0 to 5 following the invasion, and remained near zero until the blockade was lifted in September 2022.

394
chapters/C5_Object_Detection.qmd Normal file
View File

@@ -0,0 +1,394 @@
---
title: Object Detection
---
The Ship Detection tutorial explored a use case in which we might want to monitor the activity of ships in a particular location. That was a fairly straightforward task: the sea is very flat, and ships (especially large cargo and military vessels) protrude significantly. Using radar imagery, we could just set a threshold because if anything on the water is reflecting radio waves, it's probably a ship.
One shortcoming of this approach is that it doesn't tell us what *kind* of ship we've detected. Sure, you could use the shape and size to distinguish between a fishing vessel and an aircraft carrier. But what about ships of similar sizes? Or what if you wanted to use satellite imagery to identify things other than ships, like airplanes, cars, or bridges? This sort of task-- called **"object detection"** is a bit more complicated.
In this tutorial, we'll be using a deep learning model called **YOLOv5** to detect objects in satellite imagery. We'll be training the model on a custom dataset, and then using it to dynamically identify objects in satellite imagery of different resolutions pulled from Google Earth Engine. The tutorial is broken up into three sections:
1. Object detection in satellite imagery
2. Training a deep learning model on a custom dataset
3. Dynamic inference using Google Earth Engine
Unlike previous tutorials which used the GEE JavaScript API, **this one will utilize Python**; this is because these sorts of deep learning models aren't available in GEE natively yet. By the end, we'll be able to generate images such as the one below:
:::{.column-screen}
![](images/obj_det2.jpg)
:::
## Object Detection in Satellite Imagery
Object detection in satellite imagery has a variety of useful applications.
There's the needle-in-a-haystack problem of needing to monitor a large area for a small number of objects. Immediately prior to the invasion of Ukraine, for example, a number of articles emerged showing Russian military vehicles and equipment popping up in small clearings in the forest near the border with Ukraine. Many of these deployments were spotted by painstakingly combing through high resolution satellite imagery, looking for things that look like trucks. One problem with this approach is that you need to know roughly where to look. The second, and more serious problem, is that you need to be on the lookout in the first place. Object detection, applied to satellite imagery, can automatically comb through vast areas and identify objects of interest. If planes and trucks start showing up in unexpected places, you'll know about it.
Perhaps you're not monitoring that large of an area, but you want frequent updates about what's going on. What sorts of objects (planes, trucks, cars, etc.) are present? How many of each? Where are they located? Instead of having to manually look through new imagery as it becomes available, you could have a model automatically analyze new collections and output a summary.
### YOLOv5
Object detection is a fairly complicated task, and there are a number of different approaches to it. In this tutorial, we'll be using a model called **YOLOv5**. YOLO stands for **You Only Look Once**, and it's a model that was developed by [Joseph Redmon](https://pjreddie.com/) et. al., and the full paper detailing the model can be found [here](https://arxiv.org/abs/1506.02640).
The YOLOv5 model is a **convolutional neural network** (CNN), which is a type of deep learning model. CNNs are very good at identifying patterns in images, particularly in small regions of images. This is important for object detection, because we want to be able to identify objects even if they're partially obscured by other objects.
YOLO works by chopping an image up into a grid, and then predicting the location and size of objects in each grid cell:
![](images/yolo.jpg)
It learns the locations of these objects by training on a dataset of images in which each object is indicated by a **bounding box**. Then, when it's shown a new image, it will attempt to predict bounding boxes around the objects in that image. The standard YOLO model is trained on the [COCO dataset](https://cocodataset.org/#home), which contains over 200,000 images of 80 different objects ranging from people to cars to dogs. YOLO models pre-trained on this dataset work great out of the box to detect objects in videos, photographs, and live streams. But the nature of the objects we're interested in is a bit different.
Luckily, we can simply re-train the YOLOv5 model on datasets of labeled satellite imagery. The rest of this tutorial will walk through the process of training YOLOv5 on a custom dataset, and then using it to dynamically identify objects in satellite imagery pulled from Google Earth Engine.
## Training
The process of re-training the YOLOv5 model on satellite imagery is fairly straightforward and can be accomplished in just three steps; first, we're going to clone the YOLOv5 repository which contains the model code and the training scripts. Then, we'll download a dataset of satellite imagery and labels from Roboflow, and finally, we'll train the model on that dataset.
Let's start by cloning the YOLOv5 repository. Note: we'll be using a fork of the original repository that I've modified to include some pre-trained models that we'll be using later on.
```python
!git clone https://github.com/oballinger/yolov5_RS # clone repo
%cd yolov5_RS # change directory to repo
%pip install -qr requirements.txt # install dependencies
%pip install -q roboflow # install roboflow
import torch # install pytorch
import os # for os related operations
from IPython.display import Image, clear_output # to display images
```
Once we've downloaded the YOLOv5 repository, we'll need to download a dataset of labeled satellite imagery. For this example, we're going to stick with ship detection as our use case, but expand upon it. We want to be able to distinguish between different types of ships, and we want to use freely-available satellite imagery.
To that end, we'll be using [this dataset](https://universe.roboflow.com/ibl-huczk/ships-2fvbx), which contains 3400 labeled images taken from Sentinel-2 (10m/px) and PlanetScope (3m/px) satellites. Ships in these images are labeled by drawing an outline around them:
![](images/sample_training_ships.jpg)
The image above shows three ships and what is known as an STS -- a "Ship-To-Ship" transfer -- which is when a ship is transferring cargo to another ship. There are a total of seven classes of ship in this dataset:
![](images/label_freq.jpg)
This dataset can be downloaded directly from Roboflow using the following code:
```python
from roboflow import Roboflow
rf = Roboflow(api_key="<YOUR API KEY>")
project = rf.workspace('ibl-huczk').project("ships-2fvbx")
dataset = project.version("1").download("yolov5")
```
You'll need to get your own API key from Roboflow, which you can do [here](https://app.roboflow.com/account/api), and insert it in the second line of code. Roboflow is a platform for managing and training deep learning models on custom datasets. It's free to use for up to three projects, and hosts a large number of datasets that you can use to train your models. To use a different dataset, you can simply change the project name and version number in the second and third lines of code.
Finally, we can train our YOLOv5 model on the dataset we just downloaded in just one line of code:
```python
!python train.py --data {dataset.location}/data.yaml --batch 32 --cache
```
This should take about an hour.
### Accuracy Assessment
Using Tensorboard, we can log the performance of our model over the course of the training process:
:::{.column-page}
<iframe src='https://tensorboard.dev/experiment/Yiyl7AsoQcyJ3uw699CR8A/#scalars' width='100%' height='700px'></iframe>
:::
One metric in particular, **mAP 0.5**, is a good indicator of how well our model is performing. We can see it increasing rapidly at first, and then leveling off after around 30 epochs of training. The rest of this subsection will explain what exactly the mAP 0.5 value represents in this context. If you're interested in training your own model at some point, the rest of this subsection will be of interest. If you're just interested in deploying a pre-trained model, you can skip ahead to the next subsection.
In the past when we've worked on machine learning projects (for example in the makeshift refinery identifion tutorial), our training and validation data was a set of points -- geographic coordinates -- which we labeled as either being a refinery or not. Calculating the accuracy of that model was fairly straightforward, since predictions were either true positives, true negatives, false positives or false negatives.
This is slightly more complicated for object detection. We're not going pixel-by-pixel and trying to say "this is a ship" or "this is not a ship." Instead, we're looking at a larger image, and trying to draw boxes around the ships. The problem is that there are many ways to draw a box around a ship. The image below shows the labels used in our training data to indicate the location of ships.
![](images/val_batch0_labels.jpg)
![](images/val_batch0_pred.jpg)
The predicted bounding boxes are very close to the actual bounding boxes, but they're not exactly the same. The first step in evaluating the performance of our model is to determine how close the predicted boxes are to the actual boxes. We can do this by calculating the **intersection over union** (IoU) of the predicted and actual boxes. This is essentially a measure of how much overlap there is between the the predicted and actual boxes:
![Intersection over Union](images/iou.png){height=400px}
The IoU is a value between 0 and 1, where 0 means that the boxes don't overlap at all, and 1 means that the boxes overlap perfectly. Now we can set a threshold value for the IoU, and say that if the IoU is greater than that threshold, then we'll count that as a correct prediction. Now that we can classify a prediction as correct or incorrect, we can calculate two important metrics:
$$\text{Precision} = \frac{\text{True Positives}}{\text{True Positives} + \text{False Positives}}$$
This is the proportion of positive identifications that are actually correct. If my model detects 100 ships and 90 of them are actually ships, then my precision is 90%.
$$\text{Recall} = \frac{\text{True Positives}}{\text{True Positives} + \text{False Negatives}}$$
This is the proportion of actual positives that are identified correctly. If there are 100 ships in the image, and my model detects 90 of them, then my recall is 90%.
These two metrics are inversely related; I could easily get 100% recall by drawing lots of boxes everywhere to increase my chances of detecting all the ships. Conversely, I could get 100% precision by being extremely conservative and just drawing one or two boxes around the ships I'm most confident about. The key is to maximize *both*: we want our model to be sensitive enough to detect as many ships as possible (high recall), but also precise enough to only draw boxes around the ships that are actually there (high precision). Researchers find this balance using a **Precision-Recall curve** (PR curve), which plots precision on the y-axis and recall on the x-axis. Below is the Precision-Recall curve for our final model, for each class:
![Precision-Recall curve from the best ](images/pr_curve.png)
Starting from the top left corner, we set a very high confidence threshold: precision is 1, meaning that every box we draw is a ship, but recall is near 0 meaning that we're not detecting any ships. As we lower the confidence threshold, we start to detect more ships, but we also start to draw boxes around things that aren't ships. Towards the middle of the curve, we're detecting most of the ships, but we're also drawing boxes around a lot of false positives. Towards the bottom right corner, we're detecting all the ships, but we're also generating lots of false positives.
The goal is to find the point on the curve where precision and recall are both high; the closer the peak of our curve is to the top right corner, the better. A perfect model would touch the top right corner: it would have precision of 1 and recall of 1, detecting all of the ships without making any false positives. The area under this curve is called the **Average Precision** (AP), and is a measure of how close the curve is to the top right corner. The perfect model would have an AP of 1.
Some classes have a very high AP -- the value for the Aircraft Carrier class is 0.995, which is very high (though this could be down to the fact that we have a relatively small number of images with aircraft carriers in them). Ship-To-Ship (STS) transfer operations also have a high AP, at 0.951. However, other classes -- notably the "Ship" class -- have a low AP. This may be because the "Ship" class is a catch-all for any ship that doesn't fit into one of the other classes, so it encompasses lots of weird looking ships.
Finally, the **mean Average Precision** (mAP) is the average of the AP for each class, shown as the thick blue line above. Remember, all of this is premised on using a 0.5 threshold in the overlap (IoU) between our predicted boxes and the labels, which is why the final metric is called **mAP 0.5**. The mAP 0.5 for our model is 0.775, which is pretty good.
This number is very useful when training a model in several different ways using the same dataset, in order to select the best performing one. It's not that useful for comparing models trained on different datasets, since the mAP 0.5 is dependent on the number of classes in the dataset and the nature of those classes. For example, in the next section we'll be using a different model trained on the DOTA dataset which has a mAP 0.5 of around 0.68, largely due to the fact that it has around twice as many classes and many of them are similar to each other.
## Inference
Now that we've got a trained model, we can use it to conduct object detection on new images. we'll build a data processing pipeline in three steps by:
1. Loading our trained model
2. Creating an interactive map to define the area we want to analyze.
3. Defining a function to run object detection within this area.
### 1. Loading a trained model
During the training process, YOLO is iteratively tweaking the model to try to maximize mAP 0.5. It automatically saves the best version of the model in the following location: `YOLOv5_RS/runs/train/exp/weights/best.pt`. You can save this file for later use, which I have done in case you just want to use this model without having to train it yourself. I've also included several other pre-trained models which you can find in the `YOLOv5_RS/weights/` directory, including:
* `lowres_ships.pt`: the model we just trained on Sentinel-2 imagery.
* `aircraft.pt`: trained on the high resolution [Airbus Aircraft Detection Dataset](https://www.kaggle.com/datasets/airbusgeo/airbus-aircrafts-sample-dataset).
* `general.pt`: trained on the [DOTA dataset](https://captain-whu.github.io/DOTA/dataset.html) by [Kevin Guo](https://github.com/KevinMuyaoGuo/yolov5s_for_satellite_imagery#readme). This model works great on high resolution satellite imagery, and can detect the following classes: plane, ship, storage tank, baseball diamond, tennis court, basketball court, ground track field, harbor, bridge, large vehicle, small vehicle, helicopter, roundabout, soccer field, swimming pool, container crane, airport and helipad.
So far, we've trained a model to detect ships in Sentinel-2 imagery. But to show the versatility of this general approach, the rest of this tutorial will load up the `general.pt` model, and use it to detect a wide range of aircraft in high resolution imagery.
### 2. Loading the input imagery
To get started with object detection on satellite imagery using these pre-trained models, we need to define an Area of Interest (AOI) and load satellite imagery. We'll do this by accessing Google Earth Engine from the Python notebook we're working in, and creating an interactive map that will let us draw an AOI for analysis.
First, we first need to import a few packages:
```python
!pip install geemap -q
import pandas as pd
import ee
import geemap
import requests
from PIL import Image
from PIL import ImageDraw
from io import BytesIO
import torch
import PIL
```
Once we've done this, we'll also need to log in to Google Earth Engine using its Python API in order to access the satellite imagery. Running these two lines of code will generate a prompt with instructions; you have to click the link, confirm that you give the notebook permission to access your Earth Engine account, and paste the authentication code in the provided dialogue box.
```python
ee.Authenticate()
ee.Initialize()
```
Great, now we can load high resolution imagery from the National Agriculture Imagery Program (NAIP) and create an interactive map. For this example, I'm centering the map on the [Davis-Monthan Airplane Boneyard](https://en.wikipedia.org/wiki/309th_Aerospace_Maintenance_and_Regeneration_Group). This is where the US Air force retires and restores aircraft, so it will have lots of airplanes of different kinds for us to identify.
First, we want to define a function called `detect` that will accept four arguments:
* `input`: the satellite imagery we want to analyze.
* `visParams`: a dictionary of visualization parameters for the imagery.
* `weight`: the name of the pre-trained model we want to use.
* `labels`: a boolean indicating whether we want to display the labels on the processed image.
```python
def detect(input, visParams, weight, labels=True):
# Get the AOI from the map
aoi = ee.FeatureCollection(Map.draw_features)
mapScale=Map.getScale()
# Visualize the raster in Earth Engine and get a download URL
image_url=input.visualize(bands=visParams['bands'], max=visParams['max']).getThumbURL({"region":aoi.geometry(), 'scale':mapScale})
# Load the image into a PIL image
response = requests.get(image_url)
img = Image.open(BytesIO(response.content))
# Load the model
model =torch.hub.load('.','custom', path='weights/{}.pt'.format(weight),source='local',_verbose=False)
# Run inference
results = model(img)
# Count the number of detections
counts=pd.DataFrame(results.pandas().xyxy[0].groupby('name').size()).reset_index().rename(columns={0:'count','name':'detected'}).set_index('count')
# Display the results
results.show(labels=labels)
# Print the number of detections and the date of the image
print(ee.Date(input.get('system:time_start')).format("dd-MM-yyyy").getInfo())
print(counts)
return counts
```
Now, we can load the NAIP imagery and create an interactive map.
```python
# load the past 10 years of NAIP imagery
naip = ee.ImageCollection('USDA/NAIP/DOQQ').filter(ee.Filter.date('2012-01-01', '2022-01-01'))
# set some thresholds
trueColorVis = {
'bands':['R', 'G', 'B'],
'min': 0,
'max': 300,
};
# initialize our map
Map = geemap.Map()
Map.setCenter(-110.84,32.16,17)
Map.addLayer(naip.first(), trueColorVis, "naip")
Map
```
This will generate a small map with some drawing tools on the left side. We can use these tools to draw a polygon around the area we want to analyze. Use the drawing tools to draw a rectangle around an area of interest.
Finally, we can run the detection on the imagery. We'll do this by iterating through the collection of images, and running the `detect` function on each one. We'll also store the results in a dataframe so we can analyze them later.
``` python
# Get the polygon we just drew on the map
aoi=ee.FeatureCollection(Map.draw_features)
# Get a list of all the images in the collection
naip_list=naip.filterBounds(aoi).toList(naip.size())
# Iterate through the list of images and run detection on each one
for num in range(0,(img_list.size()).getInfo()):
detect(ee.Image(naip_list.get(num)), trueColorVis,'general',labels=False)
df=df.append(detection) # store the results in a dataframe
```
Below is the result of the detection on the latest image in the collection:
:::{.column-screen}
![Davis-Monthan Airplane Boneyard, Tucson AZ. 32.139498, -110.868549](images/boneyard.jpg)
:::
This image shows a remarkable degree of accuracy being achieved by our model. Inference took just 822.2 milliseconds, and it seems to be doing pretty well. The model identifies over 100 different kinds of aircraft (orange boxes) of many shapes and sizes, civilian and military, without missing a single one. It also identifies around 20 different types of helicopter (blue boxes) in the top right and even spots the cars on the highway and in the parking lots (red boxes). It's not perfect -- it thinks there's a ship in the bottom left corner near the shed (yellow box); in reality this appears to be half of a plane's fuselage, an understandable mistake given how long it took *me* to figure out what it was.
<!--
Even though we trained our model on Sentinel-2 imagery (10 meters per pixel), it can still be used on imagery from different satellites as long as they have a broadly similar resolution. A ship in PlanetScope imagery (3 meters per pixel) will look roughly similar to a ship in Sentinel-2 imagery. Using PlanetScope has another big advantage over Sentinel-2 beyond its higher spatial resolution: it has a much higher revisit rate (daily instead of 5 days). Though *downloading* PlanetScope imagery isn't free, you *can* generate a timelapse image of any area on Earth using Planet's [Planet Stories](https://www.planet.com/stories/create) tool. Simply create a free account and follow the instructions to generate a timelapse of an area of interest. You can then download the timelapse video and use it as input to our model.
Once you've done this, you can run the following line of code to automatically identify ships in the timelapse video:
{{< video images/mikolayiv.mp4 >}}
![](../images/mikolayiv.mp4)
-->

60
chapters/_quarto.yml Normal file
View File

@@ -0,0 +1,60 @@
project:
type: book
output-dir: ../docs
book:
title: "Remote Sensing \n for OSINT"
author: "[Dr. Ollie Ballinger](https://oballinger.github.io)"
date: "03/15/2023"
chapters:
- part: "A. Introduction"
chapters:
- index.qmd
- A2_Remote_Sensing.qmd
- A3_Data_Acquisition.qmd
- part: "B. Google Earth Engine"
chapters:
- B1_Getting_Started.qmd
- B2_Interpreting_Images.qmd
- B3_Image_Series.qmd
- B4_Vectors_Tables.qmd
- part: "C. Case Studies"
chapters:
- C1_Lights.qmd
- C2_Refineries.qmd
- C3_Blast.qmd
- C4_Ships.qmd
- C5_Object_Detection.qmd
repo-url: https://github.com/oballinger/RS4OSINT/
google-analytics: G-RK9ZLZQ6GL
repo-actions: [edit]
downloads: [pdf, epub]
sharing: [twitter, facebook]
favicon: ../favicon.ico
sidebar:
logo: ../images/logo_white.png
collapse-level: 1
format:
html:
theme:
dark: darkly
light: cosmo
code-copy: true
code-overflow: wrap
highlight-style: ../monokai.theme
linkcolor: "#34a832"
number-sections: false
pdf:
documentclass: scrreprt
header-includes:
- \makeatletter
- \@addtoreset{chapter}{part}
- \makeatother
epub:
cover-image: ../images/logo_black.png

233
chapters/index.log Normal file
View File

@@ -0,0 +1,233 @@
This is XeTeX, Version 3.141592653-2.6-0.999994 (TeX Live 2022) (preloaded format=xelatex 2022.10.26) 17 APR 2023 11:17
entering extended mode
restricted \write18 enabled.
%&-line parsing enabled.
**index.tex
(./index.tex
LaTeX2e <2022-06-01> patch level 5
L3 programming layer <2022-09-28> (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/scrreprt.cls
Document Class: scrreprt 2022/10/12 v3.38 KOMA-Script document class (report)
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/scrkbase.sty
Package: scrkbase 2022/10/12 v3.38 KOMA-Script package (KOMA-Script-dependent basics and keyval usage)
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/scrbase.sty
Package: scrbase 2022/10/12 v3.38 KOMA-Script package (KOMA-Script-independent basics and keyval usage)
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/scrlfile.sty
Package: scrlfile 2022/10/12 v3.38 KOMA-Script package (file load hooks)
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/scrlfile-hook.sty
Package: scrlfile-hook 2022/10/12 v3.38 KOMA-Script package (using LaTeX hooks)
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/scrlogo.sty
Package: scrlogo 2022/10/12 v3.38 KOMA-Script package (logo)
))) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/graphics/keyval.sty
Package: keyval 2022/05/29 v1.15 key=value parser (DPC)
\KV@toks@=\toks16
)
Applying: [2021/05/01] Usage of raw or classic option list on input line 252.
Already applied: [0000/00/00] Usage of raw or classic option list on input line 368.
)) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/tocbasic.sty
Package: tocbasic 2022/10/12 v3.38 KOMA-Script package (handling toc-files)
\scr@dte@tocline@numberwidth=\skip47
\scr@dte@tocline@numbox=\box51
)
Package tocbasic Info: omitting babel extension for `toc'
(tocbasic) because of feature `nobabel' available
(tocbasic) for `toc' on input line 137.
Class scrreprt Info: You've used standard option `oneside'.
(scrreprt) This is correct!
(scrreprt) Internally I'm using `twoside=false'.
(scrreprt) If you'd like to set the option with \KOMAoptions,
(scrreprt) you'd have to use `twoside=false' there
(scrreprt) instead of `oneside', too.
Class scrreprt Info: File `scrsize11pt.clo' used instead of
(scrreprt) file `scrsize11.clo' to setup font sizes on input line 2614.
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/scrsize11pt.clo
File: scrsize11pt.clo 2022/10/12 v3.38 KOMA-Script font size class option (11pt)
) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/koma-script/typearea.sty
Package: typearea 2022/10/12 v3.38 KOMA-Script package (type area)
\ta@bcor=\skip48
\ta@div=\count181
Package typearea Info: You've used standard option `letterpaper'.
(typearea) This is correct!
(typearea) Internally I'm using `paper=letter'.
(typearea) If you'd like to set the option with \KOMAoptions,
(typearea) you'd have to use `paper=letter' there
(typearea) instead of `letterpaper', too.
Package typearea Info: You've used standard option `oneside'.
(typearea) This is correct!
(typearea) Internally I'm using `twoside=false'.
(typearea) If you'd like to set the option with \KOMAoptions,
(typearea) you'd have to use `twoside=false' there
(typearea) instead of `oneside', too.
\ta@hblk=\skip49
\ta@vblk=\skip50
\ta@temp=\skip51
\footheight=\skip52
Package typearea Info: These are the values describing the layout:
(typearea) DIV = 11
(typearea) BCOR = 0.0pt
(typearea) \paperwidth = 614.295pt
(typearea) \textwidth = 446.76004pt
(typearea) DIV departure = -14%
(typearea) \evensidemargin = 11.49748pt
(typearea) \oddsidemargin = 11.49748pt
(typearea) \paperheight = 794.96999pt
(typearea) \textheight = 582.20026pt
(typearea) \topmargin = -37.40001pt
(typearea) \headheight = 17.0pt
(typearea) \headsep = 20.40001pt
(typearea) \topskip = 11.0pt
(typearea) \footskip = 47.6pt
(typearea) \baselineskip = 13.6pt
(typearea) on input line 1767.
)
\c@part=\count182
\c@chapter=\count183
\c@section=\count184
\c@subsection=\count185
\c@subsubsection=\count186
\c@paragraph=\count187
\c@subparagraph=\count188
\scr@dte@chapter@maxnumwidth=\skip53
Class scrreprt Info: using compatibility default `afterindent=bysign'
(scrreprt) for `\chapter on input line 5902.
\scr@dte@section@maxnumwidth=\skip54
Class scrreprt Info: using compatibility default `runin=bysign'
(scrreprt) for `\section on input line 5913.
Class scrreprt Info: using compatibility default `afterindent=bysign'
(scrreprt) for `\section on input line 5913.
\scr@dte@part@maxnumwidth=\skip55
Class scrreprt Info: using compatibility default `afterindent=true'
(scrreprt) for `\part on input line 5922.
\scr@dte@subsection@maxnumwidth=\skip56
Class scrreprt Info: using compatibility default `runin=bysign'
(scrreprt) for `\subsection on input line 5932.
Class scrreprt Info: using compatibility default `afterindent=bysign'
(scrreprt) for `\subsection on input line 5932.
\scr@dte@subsubsection@maxnumwidth=\skip57
Class scrreprt Info: using compatibility default `runin=bysign'
(scrreprt) for `\subsubsection on input line 5942.
Class scrreprt Info: using compatibility default `afterindent=bysign'
(scrreprt) for `\subsubsection on input line 5942.
\scr@dte@paragraph@maxnumwidth=\skip58
Class scrreprt Info: using compatibility default `runin=bysign'
(scrreprt) for `\paragraph on input line 5953.
Class scrreprt Info: using compatibility default `afterindent=bysign'
(scrreprt) for `\paragraph on input line 5953.
\scr@dte@subparagraph@maxnumwidth=\skip59
Class scrreprt Info: using compatibility default `runin=bysign'
(scrreprt) for `\subparagraph on input line 5963.
Class scrreprt Info: using compatibility default `afterindent=bysign'
(scrreprt) for `\subparagraph on input line 5963.
\abovecaptionskip=\skip60
\belowcaptionskip=\skip61
\c@pti@nb@sid@b@x=\box52
Package tocbasic Info: omitting babel extension for `lof'
(tocbasic) because of feature `nobabel' available
(tocbasic) for `lof' on input line 7140.
\scr@dte@figure@maxnumwidth=\skip62
\c@figure=\count189
Package tocbasic Info: omitting babel extension for `lot'
(tocbasic) because of feature `nobabel' available
(tocbasic) for `lot' on input line 7157.
\scr@dte@table@maxnumwidth=\skip63
\c@table=\count190
Class scrreprt Info: Redefining `\numberline' on input line 7328.
\bibindent=\dimen138
) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/amsmath/amsmath.sty
Package: amsmath 2022/04/08 v2.17n AMS math features
\@mathmargin=\skip64
For additional information on amsmath, use the `?' option.
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/amsmath/amstext.sty
Package: amstext 2021/08/26 v2.01 AMS text
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/amsmath/amsgen.sty
File: amsgen.sty 1999/11/30 v2.0 generic functions
\@emptytoks=\toks17
\ex@=\dimen139
)) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/amsmath/amsbsy.sty
Package: amsbsy 1999/11/29 v1.2d Bold Symbols
\pmbraise@=\dimen140
) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/amsmath/amsopn.sty
Package: amsopn 2022/04/08 v2.04 operator names
)
\inf@bad=\count191
LaTeX Info: Redefining \frac on input line 234.
\uproot@=\count192
\leftroot@=\count193
LaTeX Info: Redefining \overline on input line 399.
LaTeX Info: Redefining \colon on input line 410.
\classnum@=\count194
\DOTSCASE@=\count195
LaTeX Info: Redefining \ldots on input line 496.
LaTeX Info: Redefining \dots on input line 499.
LaTeX Info: Redefining \cdots on input line 620.
\Mathstrutbox@=\box53
\strutbox@=\box54
LaTeX Info: Redefining \big on input line 722.
LaTeX Info: Redefining \Big on input line 723.
LaTeX Info: Redefining \bigg on input line 724.
LaTeX Info: Redefining \Bigg on input line 725.
\big@size=\dimen141
LaTeX Font Info: Redeclaring font encoding OML on input line 743.
LaTeX Font Info: Redeclaring font encoding OMS on input line 744.
\macc@depth=\count196
LaTeX Info: Redefining \bmod on input line 905.
LaTeX Info: Redefining \pmod on input line 910.
LaTeX Info: Redefining \smash on input line 940.
LaTeX Info: Redefining \relbar on input line 970.
LaTeX Info: Redefining \Relbar on input line 971.
\c@MaxMatrixCols=\count197
\dotsspace@=\muskip16
\c@parentequation=\count198
\dspbrk@lvl=\count199
\tag@help=\toks18
\row@=\count266
\column@=\count267
\maxfields@=\count268
\andhelp@=\toks19
\eqnshift@=\dimen142
\alignsep@=\dimen143
\tagshift@=\dimen144
\tagwidth@=\dimen145
\totwidth@=\dimen146
\lineht@=\dimen147
\@envbody=\toks20
\multlinegap=\skip65
\multlinetaggap=\skip66
\mathdisplay@stack=\toks21
LaTeX Info: Redefining \[ on input line 2953.
LaTeX Info: Redefining \] on input line 2954.
) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/amsfonts/amssymb.sty
Package: amssymb 2013/01/14 v3.01 AMS font symbols
(/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/amsfonts/amsfonts.sty
Package: amsfonts 2013/01/14 v3.01 Basic AMSFonts support
\symAMSa=\mathgroup4
\symAMSb=\mathgroup5
LaTeX Font Info: Redeclaring math symbol \hbar on input line 98.
LaTeX Font Info: Overwriting math alphabet `\mathfrak' in version `bold'
(Font) U/euf/m/n --> U/euf/b/n on input line 106.
)) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/generic/iftex/iftex.sty
Package: iftex 2022/02/03 v1.0f TeX engine tests
) (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/unicode-math/unicode-math.sty (/Users/ollieballinger/Library/TinyTeX/texmf-dist/tex/latex/l3kernel/expl3.sty
Package: expl3 2023-02-07 L3 programming layer (loader)
! LaTeX Error: Mismatched LaTeX support files detected.
(LaTeX) Loading 'expl3.sty' aborted!
(LaTeX)
(LaTeX) The L3 programming layer in the LaTeX format
(LaTeX) is dated 2022-09-28, but in your TeX tree the files require
(LaTeX) at least 2023-02-07.
For immediate help type H <return>.
...
l.77 \ExplLoaderFileDate{expl3.sty}}
%
Here is how much of TeX's memory you used:
4567 strings out of 477723
94617 string characters out of 5841548
564227 words of memory out of 5000000
25658 multiletter control sequences out of 15000+600000
469267 words of font info for 29 fonts, out of 8000000 for 9000
14 hyphenation exceptions out of 8191
108i,1n,106p,10600b,270s stack positions out of 10000i,1000n,20000p,200000b,200000s
No pages of output.

35
chapters/index.qmd Normal file
View File

@@ -0,0 +1,35 @@
# Overview
The analysis of satellite imagery is a foundational element of open source investigations. Its quality, quantity and availability has increased dramatically in the past decade. Capabilities and insights that were once only available to governments are now accessible to the general public. Satellite imagery is being used to collect evidence of genocide and other war crimes in [Ukraine](https://www.nbcnews.com/science/science-news/ukraine-satellites-war-crimes-rcna26291), [Nigeria](https://www.amnesty.org/en/latest/news/2016/04/nigeria-military-cover-up-of-mass-slaughter-at-zaria-exposed/), [Burundi](https://www.amnesty.org/en/latest/news/2016/01/burundi-satellite-evidence-supports-witness-accounts-of-mass-graves/), [Cameroon](https://www.amnesty.org/en/latest/news/2021/07/cameroon-satellite-images-reveal-devastation-in-anglophone-regions-2/), [the DRC](https://www.aaas.org/resources/satellite-imagery-assessment-forced-relocations-near-luiswishi-mine), [South Sudan](https://gsp.yale.edu/case-studies/sudan/maps-satellite-images/other-darfur-satellite-imagery), [Papua](https://gsp.yale.edu/resources/maps-satellite-images/papua), and [Venezuela](https://www.hrw.org/report/2016/04/04/unchecked-power/police-and-military-raids-low-income-and-immigrant-communities). It has been used to [monitor environmental degradation](https://www.theguardian.com/environment/2016/mar/02/new-satellite-mapping-a-game-changer-against-illegal-logging) and hold extractive industries to account from [Iraq](https://www.bellingcat.com/resources/2021/04/15/what-oil-satellite-technology-and-iraq-can-tell-us-about-pollution/) to [Guatemala](https://www.planet.com/pulse/the-observatory-of-extractive-industries-oie-shines-a-light-on-the-mining-industry-using-planets-satellite-data/). The ability to analyze satellite imagery is a critical skill for anyone interested in open source investigations.
Though no-code platforms such as Sentinelhub have been invaluable in allowing the OSINT community to access and process satellite imagery, the analytical capabilities of these platforms are limited. [Google Earth Engine (GEE)](https://earthengine.google.com/#intro) is a cloud-based platform that stores petabytes of satellite imagery from a variety of sources and allows users to perform advanced analyses on Google servers for free using a browser-based interface. This textbook is designed for investigators who want to perform more sophisticated analysis using geospatial data, and assumes no prior knowledge of coding or remote sensing (satellite imagery analysis). It is organized into two parts: an introduction to remote sensing and GEE, and a series of case studies that demonstrate how to use GEE for open source investigations.
## What is Google Earth Engine?
As geospatial datasets—particularly satellite imagery collections—increase in size, researchers are increasingly relying on cloud computing platforms such as Google Earth Engine (GEE) to analyze vast quantities of data.
GEE is free and allows users to write open-source code that can be run by others in one click, thereby yielding fully reproducible results. These features have put GEE on the cutting edge of scientific research. The following plot visualizes the number of journal articles conducted using different geospatial analysis software platforms:
![](images/WoS%20Articles.png)
Despite only being released in 2015, the number of geospatial journal articles using Google Earth Engine (shown in red above) has outpaced every other major geospatial analysis software, including ArcGIS, Python, and R in just five years. GEE applications have been developed and used to present interactive geospatial data visualizations by NGOs, Universities, the United Nations, and the European Commission. By storing and running computations on Google servers, GEE is far more accessible to those who don’t have significant local computational resources; all you need is an internet connection.
## Table of Contents
A) **Introduction**
* Two introductory chapters that provide an overview of remote sensing the different types of satellite imagery available on Google Earth Engine.
* [Remote Sensing](chapters/A2_Remote_Sensing.qmd)
* [Data Acquisition](chapters/A3_Data_Acquisition.qmd)
B) **Google Earth Engine**
* Recently, a team of over 100 scientists came together to write a book called ["Cloud-Based Remote Sensing with Google Earth Engine: Fundamentals and Applications"](https://www.eefabook.org/). It's a great resource for learning about remote sensing and Earth Engine. The material in this section is a subset of the book, edited to fit the scope of this guide. If you're interested in learning more, check out the full book.
* [Getting Started](chapters/B1_Getting_Started.qmd)
* [Interpreting Images](chapters/B2_Interpreting_Images.qmd)
* [Image Series](chapters/B3_Image_Series.qmd)
* [Vectors and Tables](chapters/B4_Vectors_Tables.qmd)
C) **Case Studies**
* A series of case studies that demonstrate how to use Google Earth Engine for open source investigations. Each case study includes a brief introduction to the topic, a step-by-step guide to using Google Earth Engine to analyze satellite imagery, and a discussion of the results.
* [War at Night](chapters/C1_Lights.qmd)
* [Refinery Identification](chapters/C2_Refineries.qmd)
* [Blast Damage Assessment](chapters/C3_Blast.qmd)
* [Ship Detection](chapters/C4_Ships.qmd)
* [Object Detection](chapters/C5_Object_Detection.qmd)

19444
chapters/index.tex Normal file
View File

File diff suppressed because it is too large Load Diff