Generating cloud masks from satellite imagery - Part I

%matplotlib inline

import pandas as pd
from datetime import datetime
from IPython.display import Image
from IPython.display import HTML
import matplotlib.pyplot as plt
import numpy as np

from sklearn.preprocessing import MinMaxScaler
from datetime import datetime as dt

import arcgis
from arcgis.gis import GIS
from arcgis.learn import MLModel, prepare_tabulardata
from arcgis.raster import Raster

from fastai.vision import *

pip install sentinelhub --upgrade

pip install s2cloudless

from s2cloudless import S2PixelCloudDetector, CloudMaskRequest, get_s2_evalscript

gis = GIS("home")

# get image
s2 = gis.content.get('fd61b9e0c69c4e14bebd50a9a968348c')
sentinel = s2.layers[0]
s2

#  extent in 3857 for NYU rain clouds
NYU_cloud_extent = {
    "xmin": -8231176.77,
    "ymin": 4967559.25,
    "xmax": -8242898.16,
    "ymax": 4973524.61,
    "spatialReference": {"wkid": 3857}
}

# NYU - The respective scene having the above area is selected
selected = sentinel.filter_by(where="(Category = 1) AND (cloudcover >=0.05)",                              
                              geometry=arcgis.geometry.filters.intersects(NYU_cloud_extent))
df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear",
                   order_by_fields="AcquisitionDate").sdf

df['AcquisitionDate'] = pd.to_datetime(df['acquisitiondate'], unit='ms')
df

# The scene is selected with the least cloud cover and extracted using the selected NYC extent
NYU_clouds_scene4 = sentinel.filter_by('OBJECTID=352211')
NYU_clouds_scene4.extent = NYU_cloud_extent
NYU_clouds_scene4

from arcgis.raster.functions import apply

s2_ms = apply(NYU_clouds_scene4, None)

pd.DataFrame(s2_ms.key_properties()['BandProperties'])

# convert the scene to a numpy array
test_array = s2_ms.export_image(f='numpy_array')

test_array.shape

(450, 1200, 13)

# Change array to float data type
test_array_sr = test_array/10000
test_array_sr = test_array_sr.astype('float32')

# Initiate sentinel cloud detector 
cloud_detector = S2PixelCloudDetector(
    threshold=0.7,
    average_over=4,
    dilation_size=2,
    all_bands=True
)

%%time
test_cloud_mask = cloud_detector.get_cloud_masks(test_array_sr)

Wall time: 2.85 s

test_cloud_mask.shape

(450, 1200)

test_cloud_mask

array([[0, 0, 0, 0, ..., 0, 0, 0, 0],
       [0, 0, 0, 0, ..., 0, 0, 0, 0],
       [0, 0, 0, 0, ..., 0, 0, 0, 0],
       [0, 0, 0, 0, ..., 0, 0, 0, 0],
       ...,
       [0, 0, 0, 0, ..., 0, 0, 0, 0],
       [0, 0, 0, 0, ..., 0, 0, 0, 0],
       [0, 0, 0, 0, ..., 0, 0, 0, 0],
       [0, 0, 0, 0, ..., 0, 0, 0, 0]], dtype=int8)

plt.figure(figsize=(25,8))
plt.imshow(test_cloud_mask)
plt.title('Clouds mask by Sentinel Cloudless')
plt.axis('off')
plt.show()

#  extent in 3857 for amazon rain clouds
cloud_extent = {
    "xmin": -6501664.396,
    "ymin": -610686.285,
    "xmax": -6490039.703,
    "ymax": -600222.721,
    "spatialReference": {"wkid": 3857}
}

# The respective scene having the above area is selected
selected = sentinel.filter_by(where="(Category = 1) AND (cloudcover >=0.05)",                              
                              geometry=arcgis.geometry.filters.intersects(cloud_extent))
df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear",
                   order_by_fields="AcquisitionDate").sdf

df['AcquisitionDate'] = pd.to_datetime(df['acquisitiondate'], unit='ms')
df

# The scene is selected with the least cloud cover and extracted using the amzon extent
amazon_clouds_scene3 = sentinel.filter_by('OBJECTID=6738997')
amazon_clouds_scene3.extent = cloud_extent
amazon_clouds_scene3

s2_amazon_ms = apply(amazon_clouds_scene3, None)

pd.DataFrame(s2_amazon_ms.key_properties()['BandProperties'])

test_array = s2_amazon_ms.export_image(f='numpy_array')
test_array_sr = test_array/10000
test_array_sr = test_array_sr.astype('float32')

# Initiate sentinel cloud detector 
cloud_detector = S2PixelCloudDetector(
    threshold=0.3,
    average_over=4,
    dilation_size=2,
    all_bands=True
)

%%time
test_cloud_mask = cloud_detector.get_cloud_masks(test_array_sr)

Wall time: 2.9 s

plt.figure(figsize=(25,8))
plt.imshow(test_cloud_mask)
plt.title('Clouds mask by Sentinel Cloudless')
plt.axis('off')
plt.show()

iowa_clouds = Raster(r"https://iservices6.arcgis.com/SMX5BErCXLM7eDtY/arcgis/rest/services/small_rgb_iowa3/ImageServer", gis=gis)

iowa_clouds.extent

{
  "xmin": 659230,
  "ymin": 4642610,
  "xmax": 686720,
  "ymax": 4660770,
  "spatialReference": {
    "wkid": 32615,
    "latestWkid": 32615
  }
}

ncols = iowa_clouds.columns
nrows = iowa_clouds.rows 
iowa_clouds.read(ncols=ncols,nrows=nrows).shape

(1816, 2749, 4)

iowa_clouds.export_image(size=[2749,1816])

iowa_clouds.name

'small_rgb_iowa3'

iowa_clouds.band_names

['Band_1', 'Band_2', 'Band_3', 'Band_4']

preprocessors = [('small_rgb_iowa3', 'small_rgb_iowa3_1', 'small_rgb_iowa3_2', 'small_rgb_iowa3_3', MinMaxScaler())]

data = prepare_tabulardata(explanatory_rasters=[iowa_clouds], preprocessors=preprocessors)

data.show_batch()

model = MLModel(data, 'sklearn.cluster.KMeans', n_clusters=3, init='k-means++', random_state=43)

model.fit()

model.show_results()

# creating the local output raster path
import tempfile
local_path=tempfile.gettempdir()
output_raster_path=local_path+r"/result"+str(dt.now().microsecond)+".tif"
output_raster_path

'C:\\Users\\Supratim\\AppData\\Local\\Temp/result668383.tif'

# predicting the cloud masks using the fitted model
pred_new = model.predict(explanatory_rasters=[iowa_clouds], 
                         prediction_type='raster',
                         output_raster_path=output_raster_path)

iowa_predicted_cloud_mask = Raster(output_raster_path)
iowa_predicted_cloud_mask.export_image(size=[2749,1816])