Skip to content

Instantly share code, notes, and snippets.

@jdbcode

jdbcode/o-c-lands-from-the-perspective-of-sentinel-2.ipynb

Last active November 10, 2020 18:28
Show Gist options
  • Save jdbcode/1d21d03bc59a08ebd72ef9b3cb550eac to your computer and use it in GitHub Desktop.
Save jdbcode/1d21d03bc59a08ebd72ef9b3cb550eac to your computer and use it in GitHub Desktop.
Download ZIP
O&C Lands from the perspective of Sentinel-2
Display the source blob
Display the rendered blob
Raw
o-c-lands-from-the-perspective-of-sentinel-2.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "O&C Lands from the perspective of Sentinel-2",
"provenance": [],
"authorship_tag": "ABX9TyOJFaecRR8XkJSr9Tgbmdoz",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/jdbcode/1d21d03bc59a08ebd72ef9b3cb550eac/o-c-lands-from-the-perspective-of-sentinel-2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "OCelzsyPHmJO"
},
"source": [
"# O&C Lands from the perspective of Sentinel-2"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "rlRU5eZCGCcW"
},
"source": [
"This notebook creates a cloud-free Sentinel-2 image composite for a region in Oregon designated as Oregon and California Railroad Revested Lands (commonly known as [O&C Lands](https://en.wikipedia.org/wiki/Oregon_and_California_Railroad_Revested_Lands)). Forest harvesting patterns have left the landscape in an unnatural checkerboard pattern clearly visible from space.\n",
"\n",
"Clouds in the S2 imagery are defined by thresholding the s2cloudless collection and shadows are defined by cloud projection intersection with low-reflectance near-infrared (NIR) pixels. See this [tutorial](https://developers.google.com/earth-engine/tutorials/community/sentinel-2-s2cloudless) for more information on cloud masking with the s2cloudless collection."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "569k4qT6E6xm"
},
"source": [
"## Connect to Earth Engine servers\n",
"\n",
"Run the following cell to initialize the Earth Engine API. The output will contain instructions on how to grant this notebook access to Earth Engine using your account."
]
},
{
"cell_type": "code",
"metadata": {
"id": "6WpMP6YG91_q"
},
"source": [
"import ee\n",
"ee.Authenticate()\n",
"ee.Initialize()"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "vDWxCbuQErBf"
},
"source": [
"## Define collection filter and cloud mask parameters\n",
"\n",
"Define parameters that are used to filter the S2 image collection and determine cloud and cloud shadow identification.\n",
"\n",
"|Parameter | Type | Description |\n",
"| :-- | :-- | :-- |\n",
"| `AOI` | `ee.Geometry` | Area of interest |\n",
"| `START_DATE` | string | Image collection start date (inclusive) |\n",
"| `END_DATE` | string | Image collection end date (exclusive) |\n",
"| `CLOUD_FILTER` | integer | Maximum image cloud cover percent allowed in image collection |\n",
"| `CLD_PRB_THRESH` | integer | Cloud probability (%); values greater than are considered cloud |\n",
"| `NIR_DRK_THRESH` | float | Near-infrared reflectance; values less than are considered potential cloud shadow |\n",
"| `CLD_PRJ_DIST` | float | Maximum distance (km) to search for cloud shadows from cloud edges |\n",
"| `BUFFER` | integer | Distance (m) to dilate the edge of cloud-identified objects |"
]
},
{
"cell_type": "code",
"metadata": {
"id": "vrUkgxgA90Mi"
},
"source": [
"AOI = ee.Geometry.Rectangle(-123.5481, 43.7804, -123.2676, 43.9250)\n",
"START_DATE = '2020-07-01'\n",
"END_DATE = '2020-09-01'\n",
"CLOUD_FILTER = 60\n",
"CLD_PRB_THRESH = 40\n",
"NIR_DRK_THRESH = 0.15\n",
"CLD_PRJ_DIST = 2\n",
"BUFFER = 100"
],
"execution_count": 39,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "tcYGTOVwFNh0"
},
"source": [
"## Function to build a Sentinel-2 collection\n",
"\n",
"[Sentinel-2 surface reflectance](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR) and [Sentinel-2 cloud probability](https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_CLOUD_PROBABILITY) are two different image collections. Each collection must be filtered similarly (e.g., by date and bounds) and then the two filtered collections must be joined.\n",
"\n",
"Define a function to filter the SR and s2cloudless collections according to area of interest and date parameters, then join them on the `system:index` property. The result is a copy of the SR collection where each image has a new `'s2cloudless'` property whose value is the corresponding s2cloudless image."
]
},
{
"cell_type": "code",
"metadata": {
"id": "EsuH4LDc96ju"
},
"source": [
"def get_s2_sr_cld_col(aoi, start_date, end_date):\n",
" # Import and filter S2 SR.\n",
" s2_sr_col = (ee.ImageCollection('COPERNICUS/S2_SR')\n",
" .filterBounds(aoi)\n",
" .filterDate(start_date, end_date)\n",
" .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', CLOUD_FILTER)))\n",
"\n",
" # Import and filter s2cloudless.\n",
" s2_cloudless_col = (ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')\n",
" .filterBounds(aoi)\n",
" .filterDate(start_date, end_date))\n",
"\n",
" # Join the filtered s2cloudless collection to the SR collection by the 'system:index' property.\n",
" return ee.ImageCollection(ee.Join.saveFirst('s2cloudless').apply(**{\n",
" 'primary': s2_sr_col,\n",
" 'secondary': s2_cloudless_col,\n",
" 'condition': ee.Filter.equals(**{\n",
" 'leftField': 'system:index',\n",
" 'rightField': 'system:index'\n",
" })\n",
" }))"
],
"execution_count": 40,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "EpyDlVyOIClb"
},
"source": [
"## Function to build and apply a cloud/shadow mask"
]
},
{
"cell_type": "code",
"metadata": {
"id": "-vRZvVKs9_x8"
},
"source": [
"def apply_cld_shdw_mask(img):\n",
" # Get s2cloudless image, subset the probability band.\n",
" cld_prb = ee.Image(img.get('s2cloudless')).select('probability')\n",
"\n",
" # Condition s2cloudless by the probability threshold value.\n",
" is_cloud = cld_prb.gt(CLD_PRB_THRESH)\n",
"\n",
" # Identify water pixels from the SCL band, invert.\n",
" not_water = img.select('SCL').neq(6)\n",
"\n",
" # Identify dark NIR pixels that are not water (potential cloud shadow pixels).\n",
" SR_BAND_SCALE = 1e4\n",
" dark_pixels = img.select('B8').lt(NIR_DRK_THRESH*SR_BAND_SCALE).multiply(not_water)\n",
"\n",
" # Determine the direction to project cloud shadow from clouds (assumes UTM projection).\n",
" shadow_azimuth = ee.Number(90).subtract(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE')));\n",
"\n",
" # Project shadows from clouds for the distance specified by the CLD_PRJ_DIST input.\n",
" cld_proj = (is_cloud.directionalDistanceTransform(shadow_azimuth, CLD_PRJ_DIST*10)\n",
" .reproject(**{'crs': img.select(0).projection(), 'scale': 100})\n",
" .select('distance')\n",
" .mask())\n",
"\n",
" # Identify the intersection of dark pixels with cloud shadow projection.\n",
" is_shadow = cld_proj.multiply(dark_pixels)\n",
"\n",
" # Combine cloud and shadow mask, set cloud and shadow as value 1, else 0.\n",
" is_cld_shdw = is_cloud.add(is_shadow).gt(0)\n",
"\n",
" # Remove small cloud-shadow patches and dilate remaining pixels by BUFFER input.\n",
" # 20 m scale is for speed, and assumes clouds don't require 10 m precision.\n",
" is_cld_shdw = (is_cld_shdw.focal_min(2).focal_max(BUFFER*2/20)\n",
" .reproject(**{'crs': img.select([0]).projection(), 'scale': 20})\n",
" .rename('cloudmask'))\n",
"\n",
" # Subset reflectance bands and update their masks, return the result.\n",
" return img.select('B.*').updateMask(is_cld_shdw.Not())"
],
"execution_count": 41,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "jxbA458EIy48"
},
"source": [
"## Make a cloud-free image composite and display it."
]
},
{
"cell_type": "code",
"metadata": {
"id": "BwcxGN1J-Jh0"
},
"source": [
"s2_sr_cld_col = get_s2_sr_cld_col(AOI, START_DATE, END_DATE)\n",
"s2_sr_median = s2_sr_cld_col.map(apply_cld_shdw_mask).median()"
],
"execution_count": 42,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "v3AwWOTIALpS"
},
"source": [
"from IPython.display import Image\n",
"\n",
"print('Please wait patiently, processing may take a minute...')\n",
"vis_params = {\n",
" 'bands': ['B11', 'B8', 'B4'],\n",
" 'min': 100, 'max': 4000, 'gamma': 0.7,\n",
" 'dimensions': 1024, 'region': AOI, 'crs': 'EPSG:3857',\n",
" 'format': 'PNG'\n",
"}\n",
"\n",
"Image(url = (s2_sr_median.getThumbURL(vis_params)))"
],
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment