This guide is a follow up to the UDM2 Cloud Detection notebook. Please refer to that notebook for further details on specifications and usage of the UDM2 asset. In this notebook, we apply cloud detection to a specific area of interest (AOI).
The udm2 asset is available for all PSScene4Band and PSOrthoTile items created after 2018-08-01. Therefore, our search should be limited to these items and this date range. In this notebook, we focus on PSOrthoTile imagery taken within the month of April, 2019.
In [1]:
import datetime
import json
import os
import time
from planet import api
from planet.api import filters
import rasterio
from rasterio.plot import show
One of the benefits of accurate and automated cloud detection is that it allows users to filter out images that don't meet a certain quality threshold. Planet's Data API allows users to search based on the value of the imagery metadata.
Planet's cloud detection algorithm classifies every pixel into one of six different categories, each of which has a corresponding metadata field that reflects the percentage of data that falls into the category.
| Class | Metadata field |
|---|---|
| clear | clear_percent |
| snow | snow_ice_percent |
| shadow | shadow_percent |
| light haze | light_haze_percent |
| heavy haze | heavy_haze_percent |
| cloud | cloud_percent |
The UDM2 Cloud Detection notebook provides examples for how to perform searches using the Planet Python Client command-line interface.
The following example will show how use Planet's Python Client to perform a search for PSOrthoTiles that are 70-95% clear (not totally clear so the udm has some interesting content), taken within the month of April, 2019, and within an AOI.
In [2]:
# define the date range for imagery
start_date = datetime.datetime(year=2019,month=4,day=1)
stop_date = datetime.datetime(year=2019,month=5,day=1)
# define the aoi for imagery
aoi = {"geometry": {
"type":"Polygon",
"coordinates":
[[
[25.42429478260258,1.0255377823058893],
[25.592960813580472,1.0255377823058893],
[25.592960813580472,1.1196578801254304],
[25.42429478260258,1.1196578801254304],
[25.42429478260258,1.0255377823058893]
]]}}
item_type = 'PSOrthoTile'
In [19]:
# create an api request from the search specifications
def build_ps_request(item_type, aoi_geom, start_date, stop_date):
query = filters.and_filter(
filters.geom_filter(aoi_geom),
filters.range_filter('clear_percent', gte=70),
filters.range_filter('clear_percent', lte=95),
filters.date_range('acquired', gt=start_date),
filters.date_range('acquired', lt=stop_date)
)
return filters.build_search_request(query, [item_type])
request = build_ps_request(item_type, aoi['geometry'], start_date, stop_date)
print(request)
In [4]:
# search the API for imagery
def get_api_key():
return os.environ['PL_API_KEY']
# quick check that key is defined
assert get_api_key(), "PL_API_KEY not defined."
def create_client():
return api.ClientV1(api_key=get_api_key())
def search_pl_api(request, limit=500):
client = create_client()
result = client.quick_search(request)
# note that this returns a generator
return result.items_iter(limit=limit)
items = list(search_pl_api(request))
print(len(items))
# uncomment below to see entire metadata for a PS orthotile
# print(json.dumps(items[0], indent=4))
In [5]:
item_id = items[0]['id']
item_id
Out[5]:
In addition to metadata for filtering, the udm2 asset provides a pixel-by-pixel map that identifies the classification of each pixel. See the UDM2 Cloud Detection notebook for an example map.
The udm2 structure is to use a separate band for each classification type. Band 2, for example, indicates that a pixel is snowy when its value is 1, band 3 indicates shadow and so on. "
The following Python will download the data above and then display pixels that fall into a certain classifications.
In [7]:
# activate assets
client = create_client()
assets = client.get_assets_by_id(item_type, item_id).get()
client.activate(assets["analytic"])
client.activate(assets["udm2"])
# wait until activation completes
while True:
assets = client.get_assets_by_id(item_type, item_id).get()
if "location" in assets["analytic"] and "location" in assets["udm2"]:
print('assets activated')
break
print('waiting')
time.sleep(15)
In [8]:
# start downloads
data_dir = 'data'
r1 = client.download(assets["analytic"], callback=api.write_to_file(data_dir))
r2 = client.download(assets["udm2"], callback=api.write_to_file(data_dir))
In [9]:
# wait until downloads complete
r1.wait()
r2.wait()
img_file = os.path.join(data_dir, r1.get_body().name)
udm_file = os.path.join(data_dir, r2.get_body().name)
print("image: {}".format(img_file))
print("udm2: {}".format(udm_file))
In [20]:
with rasterio.open(udm_file) as src:
shadow_mask = src.read(3).astype(bool)
cloud_mask = src.read(6).astype(bool)
show(shadow_mask, title="shadow", cmap="binary")
show(cloud_mask, title="cloud", cmap="binary")
Out[20]:
In [18]:
mask = shadow_mask + cloud_mask
show(mask, title="mask", cmap="binary")
Out[18]:
In [13]:
with rasterio.open(img_file) as src:
profile = src.profile
img_data = src.read([3, 2, 1], masked=True) / 10000.0 # apply RGB ordering and scale down
In [14]:
show(img_data, title=item_id)
Out[14]:
In [15]:
img_data.mask = mask
img_data = img_data.filled(fill_value=0)
In [16]:
show(img_data, title="masked image")
Out[16]:
The image stored in img_data now has cloudy / cloud-shadowy pixels masked out and can be saved or used for analysis.