Skip to content

Commit fa10dc4

Browse files
thuydotmthuydotm
authored
added drop_clouds, cloud-free mosaic from sentinel2 data example (#255)
1 parent 16b2b84 commit fa10dc4

File tree

3 files changed

+633
-13
lines changed

3 files changed

+633
-13
lines changed

‎examples/cloudless-mosaic-sentinel2.ipynb

+263
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,263 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Mosaic from a single multitemporal dataset\n",
8+
"\n",
9+
"\n",
10+
"The goal of this notebook is to provide an example of how to create a cloud-free mosaic from Sentinel-2 imagery over a specific area over a time period. We first use `satsearch` to search for Sentinel-2 data then combine them together using `stackstac`. A median operation will be applied to merge the images into a single layer that could be save off into Azure blob storage as COGs for later use.\n",
11+
"\n",
12+
"\n",
13+
"## 1. Sentinel-2 Dataset\n",
14+
"\n",
15+
"Satellite images (also Earth observation imagery, spaceborne photography, or simply satellite photo) are images of Earth collected by imaging satellites operated by governments and businesses around the world (see https://en.wikipedia.org/wiki/Satellite_imagery). Its major applications include Earth observation and land cover monitoring. \n",
16+
"\n",
17+
"\n",
18+
"SENTINEL-2 (https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/overview) is a wide-swath, high-resolution, multi-spectral imaging mission, supporting Copernicus Land Monitoring studies, including the monitoring of vegetation, soil and water cover, as well as observation of inland waterways and coastal areas.\n",
19+
"\n",
20+
"## 2. Imports"
21+
]
22+
},
23+
{
24+
"cell_type": "code",
25+
"execution_count": null,
26+
"metadata": {},
27+
"outputs": [],
28+
"source": [
29+
"import stackstac\n",
30+
"from satsearch import Search\n",
31+
"\n",
32+
"import xrspatial.multispectral as ms\n",
33+
"\n",
34+
"import matplotlib.pyplot as plt"
35+
]
36+
},
37+
{
38+
"cell_type": "markdown",
39+
"metadata": {},
40+
"source": [
41+
"## 3. Load Sentinel 2 data\n",
42+
"\n",
43+
"In this example, we use data from `sentinel-s2-l2a-cogs` collection within a bounding box of `[-93.112301, 29.649001, -92.075965, 30.719868]`, and the time range considered is from `2019-07-01` to `2020-06-30`. And the collected data has less than 25% cloud coverage."
44+
]
45+
},
46+
{
47+
"cell_type": "code",
48+
"execution_count": null,
49+
"metadata": {},
50+
"outputs": [],
51+
"source": [
52+
"items = Search(\n",
53+
" url=\"https://earth-search.aws.element84.com/v0\",\n",
54+
" bbox=[-93.112301, 29.649001, -92.075965, 30.719868],\n",
55+
" collections=[\"sentinel-s2-l2a-cogs\"],\n",
56+
" query={'eo:cloud_cover': {'lt': 25}},\n",
57+
" datetime=\"2019-07-01/2020-06-30\"\n",
58+
").items()\n",
59+
"\n",
60+
"len(items)"
61+
]
62+
},
63+
{
64+
"cell_type": "markdown",
65+
"metadata": {},
66+
"source": [
67+
"Let's combine all the above STAC items into a lazy xarray with following settings:\n",
68+
"- projection: epsg=32613\n",
69+
"- resolution: 100m\n",
70+
"- bands: green (B02), red (B03), blue (B04), NIR (B08), SWIR1 (B11)"
71+
]
72+
},
73+
{
74+
"cell_type": "code",
75+
"execution_count": null,
76+
"metadata": {},
77+
"outputs": [],
78+
"source": [
79+
"bands = ['B02', 'B03', 'B04', 'B08', 'B11']\n",
80+
"\n",
81+
"stack_ds = stackstac.stack(\n",
82+
" items, epsg=32613, resolution=100, assets=bands\n",
83+
")\n",
84+
"\n",
85+
"stack_ds"
86+
]
87+
},
88+
{
89+
"cell_type": "markdown",
90+
"metadata": {},
91+
"source": [
92+
"We can get a median composite for each month in the considered period of time:"
93+
]
94+
},
95+
{
96+
"cell_type": "code",
97+
"execution_count": null,
98+
"metadata": {},
99+
"outputs": [],
100+
"source": [
101+
"monthly = stack_ds.resample(time=\"MS\").median(\"time\", keep_attrs=True)\n",
102+
"monthly"
103+
]
104+
},
105+
{
106+
"cell_type": "code",
107+
"execution_count": null,
108+
"metadata": {},
109+
"outputs": [],
110+
"source": [
111+
"import dask.diagnostics\n",
112+
"with dask.diagnostics.ProgressBar():\n",
113+
" monthly = monthly.compute()"
114+
]
115+
},
116+
{
117+
"cell_type": "markdown",
118+
"metadata": {},
119+
"source": [
120+
"## 4. Cloud-free scene using median operator\n",
121+
"\n",
122+
"In this step, we use a median operation to merge all monthly images into 1 single cloud-free layer. With an assumption that, along a multitemporal stack, clouds would not persist at the same geographical position from time to time (i.e image to image), the more data we have, the higher chance of dropping clouds."
123+
]
124+
},
125+
{
126+
"cell_type": "code",
127+
"execution_count": null,
128+
"metadata": {},
129+
"outputs": [],
130+
"source": [
131+
"median_scene = monthly.median(dim=['time'])"
132+
]
133+
},
134+
{
135+
"cell_type": "markdown",
136+
"metadata": {},
137+
"source": [
138+
"With 3 bands: red, green, blue, let's see the true color using the `true_color` function from `xrspatial.multispectral module` for each separate month and the median layer."
139+
]
140+
},
141+
{
142+
"cell_type": "code",
143+
"execution_count": null,
144+
"metadata": {},
145+
"outputs": [],
146+
"source": [
147+
"bands_mapping = {v: i for i, v in enumerate(bands)}\n",
148+
"\n",
149+
"band_blue = bands_mapping['B02']\n",
150+
"band_green = bands_mapping['B03']\n",
151+
"band_red = bands_mapping['B04']"
152+
]
153+
},
154+
{
155+
"cell_type": "code",
156+
"execution_count": null,
157+
"metadata": {},
158+
"outputs": [],
159+
"source": [
160+
"months = 12\n",
161+
"imgs = []\n",
162+
"for month in range(months):\n",
163+
" # True color\n",
164+
" r = monthly[month][band_red]\n",
165+
" g = monthly[month][band_green]\n",
166+
" b = monthly[month][band_blue]\n",
167+
" img = ms.true_color(r, g, b)\n",
168+
" imgs.append(img)"
169+
]
170+
},
171+
{
172+
"cell_type": "code",
173+
"execution_count": null,
174+
"metadata": {},
175+
"outputs": [],
176+
"source": [
177+
"# Utility function for displaying images\n",
178+
"\n",
179+
"def display_images(images, columns=2, width=50, height=50):\n",
180+
" height = max(height, int(len(images)/columns) * height)\n",
181+
" plt.figure(figsize=(width, height))\n",
182+
" for i, image in enumerate(images):\n",
183+
" plt.subplot(len(images) / columns + 1, columns, i + 1)\n",
184+
" plt.imshow(image)"
185+
]
186+
},
187+
{
188+
"cell_type": "markdown",
189+
"metadata": {},
190+
"source": [
191+
"#### Monthly data"
192+
]
193+
},
194+
{
195+
"cell_type": "code",
196+
"execution_count": null,
197+
"metadata": {},
198+
"outputs": [],
199+
"source": [
200+
"# takes some time to run\n",
201+
"display_images(imgs)"
202+
]
203+
},
204+
{
205+
"cell_type": "markdown",
206+
"metadata": {},
207+
"source": [
208+
"#### Median layer"
209+
]
210+
},
211+
{
212+
"cell_type": "code",
213+
"execution_count": null,
214+
"metadata": {},
215+
"outputs": [],
216+
"source": [
217+
"median_scene = monthly.median(dim=['time'])\n",
218+
"\n",
219+
"median_red_agg = median_scene[band_red]\n",
220+
"median_green_agg = median_scene[band_green]\n",
221+
"median_blue_agg = median_scene[band_blue]\n",
222+
"\n",
223+
"median_img = ms.true_color(median_red_agg, median_green_agg, median_blue_agg)\n",
224+
"\n",
225+
"median_img"
226+
]
227+
},
228+
{
229+
"cell_type": "markdown",
230+
"metadata": {},
231+
"source": [
232+
"### References\n",
233+
"\n",
234+
"- https://stackstac.readthedocs.io/en/latest/basic.html"
235+
]
236+
}
237+
],
238+
"metadata": {
239+
"kernelspec": {
240+
"display_name": "Python 3",
241+
"language": "python",
242+
"name": "python3"
243+
},
244+
"language_info": {
245+
"codemirror_mode": {
246+
"name": "ipython",
247+
"version": 3
248+
},
249+
"file_extension": ".py",
250+
"mimetype": "text/x-python",
251+
"name": "python",
252+
"nbconvert_exporter": "python",
253+
"pygments_lexer": "ipython3",
254+
"version": "3.8.5"
255+
},
256+
"widgets": {
257+
"state": {},
258+
"version": "1.1.2"
259+
}
260+
},
261+
"nbformat": 4,
262+
"nbformat_minor": 4
263+
}

0 commit comments

Comments
 (0)