Skip to content

Instantly share code, notes, and snippets.

@jdbcode

jdbcode/irrigated-fields-sentinel-2.ipynb

Last active January 10, 2023 00:54
Show Gist options
  • Save jdbcode/fc4a2d15bf26b4635860d3296cb78411 to your computer and use it in GitHub Desktop.
Save jdbcode/fc4a2d15bf26b4635860d3296cb78411 to your computer and use it in GitHub Desktop.
Download ZIP
Irrigated fields Sentinel-2
Display the source blob
Display the rendered blob
Raw
irrigated-fields-sentinel-2.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "Irrigated fields Sentinel-2",
"provenance": [],
"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/fc4a2d15bf26b4635860d3296cb78411/irrigated-fields-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": [
"# Irrigated fields Sentinel-2"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "rlRU5eZCGCcW"
},
"source": [
"This notebook creates a cloud-free Sentinel-2 image composite for a region with irrigated fields along the Columbia river in eastern Washington.\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(-120.0433, 46.5950, -119.5036, 46.8247)\n",
"START_DATE = '2020-07-01'\n",
"END_DATE = '2020-09-01'\n",
"CLOUD_FILTER = 20\n",
"CLD_PRB_THRESH = 70 # 40\n",
"NIR_DRK_THRESH = 0.15 # 0.15\n",
"CLD_PRJ_DIST = 2\n",
"BUFFER = 75"
],
"execution_count": null,
"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_HARMONIZED) 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_HARMONIZED')\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": null,
"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": null,
"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": null,
"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': 4500, 'gamma': 0.7,\n",
" 'dimensions': 1024, 'region': AOI, 'crs': 'EPSG:3857',\n",
" 'format': 'PNG'\n",
"}\n",
"\n",
"url = s2_sr_median.getThumbURL(vis_params)\n",
"Image(url=url)"
],
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment