S2_SR_8bit_quarterly_mosaics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ea4644ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Create a truecolor or falsecolor, cloud-masked 8-bit RGB mosaic from Sentinel-2 Surface Reflectance Data.\n",
    "https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR_HARMONIZED#description\n",
    "Relies on the s2cloudless dataset: https://developers.google.com/earth-engine/tutorials/community/sentinel-2-s2cloudless\n",
    "\n",
    "Please note, for all of BC the mosaic will be approximately 64 GB. \n",
    "\n",
    "Quarters:\n",
    "q1: Jan 1 - Mar 30\n",
    "q2: Apr 1 - Jun 30\n",
    "q3: Jul 1 - Sept 30\n",
    "q4: Oct 1 - Dec 31\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "#Parameters\n",
    "area_name = 'test' # 'bc' or anything else\n",
    "\n",
    "if area_name == 'bc':\n",
    "    p = r\"C:\\Data\\Datasets\\WHSE_ADMIN_BOUNDARIES_ADM_NR_AREAS_SP_dissolved.shp\"\n",
    "else:\n",
    "    p = r\"C:\\Data\\Datasets\\test.shp\"\n",
    "\n",
    "\n",
    "q = 'q2' #q1, q2, q3, q4 or custom\n",
    "comp = 'tc' #tc (truecolor) or fc (falsecolor)\n",
    "year = '2023'\n",
    "\n",
    "##for custom define dates \"yyyy-mm-dd\"\n",
    "if q == 'custom':\n",
    "    t1 = '2023-06-15'\n",
    "    t2 = '2023-07-04'\n",
    "\n",
    "# Output parameters\n",
    "scale = 10 #in meters\n",
    "proj = 'EPSG:3005' #epsg code\n",
    "outfolder = 'mosaic_' + year + '_' + q #output google drive folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "a9d39a78",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import geemap\n",
    "import os\n",
    "import math\n",
    "from osgeo import gdal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "6aee3247",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getfiles(d,ext):\n",
    "    paths = []\n",
    "    for file in os.listdir(d):\n",
    "        if file.endswith(ext):\n",
    "            paths.append(os.path.join(d, file))\n",
    "    return(paths) \n",
    "\n",
    "###Cloud masking functions for s2cloudless \n",
    "###(https://developers.google.com/earth-engine/tutorials/community/sentinel-2-s2cloudless)\n",
    "\n",
    "# Get surface reflectance data and corresponding s2cloudless dataset, merge\n",
    "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",
    "    }))\n",
    "\n",
    "# Get cloud bands\n",
    "# Cloud probability is defined in the main body of the script, using 50%, default is 30% but I found that too aggressive\n",
    "def add_cloud_bands(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).rename('clouds')\n",
    "\n",
    "    # Add the cloud probability layer and cloud mask as image bands.\n",
    "    return img.addBands(ee.Image([cld_prb, is_cloud]))\n",
    "\n",
    "# Create cloud shadow bands\n",
    "# NIR_DRK_THRESH is defined in the main body of the script\n",
    "def add_shadow_bands(img):\n",
    "    # Identify water pixels from the SCL band.\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).rename('dark_pixels')\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 = (img.select('clouds').directionalDistanceTransform(shadow_azimuth, CLD_PRJ_DIST*10)\n",
    "        .reproject(**{'crs': img.select(0).projection(), 'scale': 100})\n",
    "        .select('distance')\n",
    "        .mask()\n",
    "        .rename('cloud_transform'))\n",
    "\n",
    "    # Identify the intersection of dark pixels with cloud shadow projection.\n",
    "    shadows = cld_proj.multiply(dark_pixels).rename('shadows')\n",
    "\n",
    "    # Add dark pixels, cloud projection, and identified shadows as image bands.\n",
    "    return img.addBands(ee.Image([dark_pixels, cld_proj, shadows]))\n",
    "\n",
    "# Merge cloud and cloud shadow mask, clean up\n",
    "# BUFFER is defined in the main body of the script, default 50m, I use 10m. Again, too aggressive\n",
    "# There is confusion between snow and cloud\n",
    "\n",
    "def add_cld_shdw_mask(img):\n",
    "    # Add cloud component bands.\n",
    "    img_cloud = add_cloud_bands(img)\n",
    "\n",
    "    # Add cloud shadow component bands.\n",
    "    img_cloud_shadow = add_shadow_bands(img_cloud)\n",
    "\n",
    "    # Combine cloud and shadow mask, set cloud and shadow as value 1, else 0.\n",
    "    is_cld_shdw = img_cloud_shadow.select('clouds').add(img_cloud_shadow.select('shadows')).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.focalMin(2).focalMax(BUFFER*2/20)\n",
    "        .reproject(**{'crs': img.select([0]).projection(), 'scale': 20})\n",
    "        .rename('cloudmask'))\n",
    "\n",
    "    # Add the final cloud-shadow mask to the image.\n",
    "    return img_cloud_shadow.addBands(is_cld_shdw)\n",
    "\n",
    "# Apply all masks\n",
    "def apply_cld_shdw_mask(img):\n",
    "    # Subset the cloudmask band and invert it so clouds/shadow are 0, else 1.\n",
    "    not_cld_shdw = img.select('cloudmask').Not()\n",
    "\n",
    "    # Subset reflectance bands and update their masks, return the result.\n",
    "    return img.select('B.*').updateMask(not_cld_shdw)\n",
    "\n",
    "## Other functions\n",
    "# Get min and max of an image, for debug \n",
    "def get_minmax(image,poly):\n",
    "    min_max = ee.Image(image).reduceRegion(**{\n",
    "        'reducer': ee.Reducer.minMax(),\n",
    "        'geometry': poly.geometry(),\n",
    "        'scale': 100,\n",
    "        'maxPixels': 1e10})\n",
    "    print(min_max.getInfo())\n",
    "\n",
    "# Convert to 8-bit\n",
    "# Apply pre-defined gamma values\n",
    "# Retain 0 as no-data value so stretch between 1 and 255\n",
    "def convert(image=None,gamma=None,maxValue=None,rgbType=None):\n",
    "    if rgbType == 'tc':\n",
    "        img = image.select(['B4','B3','B2'])\n",
    "    else:\n",
    "        img = image.select(['B11', 'B8', 'B4'])\n",
    "    \n",
    "    #reclass values less than 0 to 0\n",
    "    img1 = img.where(img.gt(maxValue),maxValue) \n",
    "    \n",
    "    #reclass values greater than 0.4 to 0.4\n",
    "    img2 = img1.where(img1.lt(0),0)\n",
    "    \n",
    "    #apply gamma function\n",
    "    eqtn = 'pow(x, 1.0/' + str(gamma) + ')'\n",
    "    gamma_img = img2.expression(eqtn,{'x':img2})\n",
    "    \n",
    "    #get max gamma value for 8-bit transformation\n",
    "    max_gamma = math.pow(maxValue,1.0/gamma)\n",
    "    \n",
    "    #Convert to 1 - 255\n",
    "    eqtn = '((x/' + str(max_gamma) + ')' +'*254)+1'\n",
    "    img_out = gamma_img.expression(eqtn,{'x':gamma_img}).int()\n",
    "    return(img_out)\n",
    "\n",
    "def add_alpha(img):\n",
    "    eqtn = \"b('vis-red') > 0 ? 255 : 0\"\n",
    "    alpha = img.expression(eqtn,{'x':img}).rename('alpha').toUint8()\n",
    "    return(img.addBands(alpha))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "36960771",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Intialize gee\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "474af0bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import vector data\n",
    "poly = geemap.shp_to_ee(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "9f4522dc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2023-04-01\n",
      "2023-06-30\n"
     ]
    }
   ],
   "source": [
    "#Look for Sentinel-2 data, apply cloud masks, create a median mosaic, clip to extent of the polygon provided \n",
    "#Quarterly mosaics: Q1 (Jan1 - Mar31), Q2 (Apr1 - Jun30), Q3 (Jul1 - Sept30), Q4 (Oct1 - Dec31)\n",
    "\n",
    "if q == 'q1':\n",
    "    START_DATE= year + '-01-01'\n",
    "    END_DATE = year + '-03-31'\n",
    "    stretch = 'winter'\n",
    "    interval = year + q\n",
    "elif q == 'q2':\n",
    "    START_DATE= year + '-04-01'\n",
    "    END_DATE = year + '-06-30'\n",
    "    stretch = 'winter'\n",
    "    interval = year + q\n",
    "elif q == 'q3':\n",
    "    START_DATE= year + '-07-01'\n",
    "    END_DATE = year + '-09-30'\n",
    "    stretch = 'summer'\n",
    "    interval = year + q\n",
    "elif q == 'q4':\n",
    "    START_DATE= year + '-10-01'\n",
    "    END_DATE = year + '-12-31'\n",
    "    stretch = 'winter'\n",
    "    interval = year + q\n",
    "elif q == 'custom':\n",
    "    START_DATE = t1\n",
    "    END_DATE = t2\n",
    "    stretch = 'summer'\n",
    "    interval = START_DATE + '_' + END_DATE\n",
    "else:\n",
    "    print('no quarter specified')\n",
    "\n",
    "## Define stretches\n",
    "if stretch == 'summer':\n",
    "    if comp == 'tc':\n",
    "        comp_str = '_b432'\n",
    "        gamma=1.5\n",
    "        maxValue=0.4\n",
    "    elif comp == 'fc':\n",
    "        comp_str = '_b1184'\n",
    "        gamma=1.5\n",
    "        maxValue=0.4\n",
    "elif stretch == 'winter':\n",
    "    if comp == 'tc':\n",
    "        comp_str = '_b432'\n",
    "        gamma=2\n",
    "        maxValue=1.0\n",
    "        rgbType=comp\n",
    "    elif comp == '_fc':\n",
    "        comp_str = '_b1184'\n",
    "        gamma=2\n",
    "        maxValue=1.0\n",
    "    \n",
    "CLOUD_FILTER = 80\n",
    "CLD_PRB_THRESH = 50\n",
    "NIR_DRK_THRESH = 0.15 \n",
    "CLD_PRJ_DIST = 1\n",
    "BUFFER = 10 \n",
    "\n",
    "viz = {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max':255}\n",
    "\n",
    "#Search Sentinel-2 Surface Reflectance collection\n",
    "s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED').filterDate(START_DATE,END_DATE).filterBounds(poly).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',CLOUD_FILTER))\n",
    "\n",
    "## Create a cloud masked mosaic\n",
    "col = get_s2_sr_cld_col(poly, START_DATE, END_DATE)\n",
    "col_wmsks = col.map(add_cld_shdw_mask)\n",
    "s2_noclouds = col_wmsks.map(apply_cld_shdw_mask).median().multiply(0.0001)\n",
    "s2_noclouds_med_vis = convert(image=s2_noclouds,gamma=gamma,maxValue=maxValue,rgbType=comp).clip(poly).visualize(**viz)\n",
    "\n",
    "#get_minmax(s2_noclouds_med_vis,poly) #for debug, very slow for large images\n",
    "\n",
    "## Add alpha band B1 > 0 == 255, name alpha\n",
    "test = add_alpha(s2_noclouds_med_vis)\n",
    "test\n",
    "\n",
    "print(START_DATE)\n",
    "print(END_DATE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4481b10d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "22f8bd10294d4557896ed6401974e3f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[54.35439399296694, -124.80707076512125], controls=(WidgetControl(options=['position', 'transparent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Display\n",
    "## Zoom in, sometimes won't show the whole provincial mosaic when zoomed out.\n",
    "## Works well with smaller areas\n",
    "\n",
    "Map = geemap.Map()\n",
    "\n",
    "#add to Map\n",
    "s2_vis = {\n",
    "  'min': 0,\n",
    "  'max': 255,\n",
    "  'bands':['vis-red', 'vis-green', 'vis-blue'],\n",
    "};\n",
    "\n",
    "Map.addLayer(s2_noclouds_med_vis, s2_vis, 'mosaic')\n",
    "Map.centerObject(poly, 10)\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0cf85468",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "s2_bc_2023q2_bcalb_10m_b432\n"
     ]
    }
   ],
   "source": [
    "##Export to google drive\n",
    "region = poly.geometry()\n",
    "outname = 's2_' + area_name + '_' + interval + '_bcalb_' + str(scale) +'m_b432'\n",
    "print(outname)\n",
    "\n",
    "task_config = {\n",
    "    'description': outname,\n",
    "    'scale': scale,  \n",
    "    'region': region,\n",
    "    'crs': proj,\n",
    "    'folder':outfolder,\n",
    "    'shardSize':256,\n",
    "    'fileDimensions':1024,\n",
    "    'maxPixels':10000000000000,\n",
    "    'skipEmptyTiles':True}\n",
    "\n",
    "task = ee.batch.Export.image.toDrive(s2_noclouds_med_vis.select(['vis-red', 'vis-green', 'vis-blue']), **task_config)\n",
    "task.start()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9193c06a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x000002B1D875E460> >"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Mosaic tiles once downloaded\n",
    "import os\n",
    "from osgeo import gdal\n",
    "\n",
    "f = r'C:\\Data\\2023\\q2\\tiles'\n",
    "outfolder = r'C:\\Data\\2023\\q2\\mosaic'\n",
    "base = 's2_bc_2023q2_bcalb_10m_b432'\n",
    "\n",
    "files = getfiles(f,'.tif')\n",
    "outfilename = os.path.join(outfolder,base+'.tif')\n",
    "gdal.Warp(outfilename,files)\n",
    "\n",
    "#set no data value to which makes it transparent \n",
    "outdata = gdal.Open(os.path.join(outfolder,base+'.tif'))\n",
    "outdata.GetRasterBand(1).SetNoDataValue(0) \n",
    "outdata = None  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ec9da44",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}