Skip to content

Instantly share code, notes, and snippets.

@RutgerK

RutgerK/slstr_cloudmask.ipynb

Created July 16, 2018 09:46
Show Gist options
  • Save RutgerK/3663fea47c03013b11a69512a3971419 to your computer and use it in GitHub Desktop.
Save RutgerK/3663fea47c03013b11a69512a3971419 to your computer and use it in GitHub Desktop.
Download ZIP
SLSTR_cloudmask
Display the source blob
Display the rendered blob
Raw
slstr_cloudmask.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"D:\\Algemeen\\VSCode\\satpy\n"
]
}
],
"source": [
"%cd D:\\Algemeen\\VSCode\\satpy"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from satpy import Scene\n",
"from satpy import find_files_and_readers\n",
"from pyresample.geometry import AreaDefinition\n",
"\n",
"import xarray as xr\n",
"import numpy as np\n",
"\n",
"import os\n",
"import datetime"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Settings"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"sensor = 'slstr'\n",
"reader = 'nc_slstr'\n",
"\n",
"base_dir = r'D:\\Algemeen\\Sentinel3\\data'\n",
"\n",
"area_slstr = AreaDefinition('tmp', 'tmp', 'tmp', dict(init='EPSG:4326'), 2273, 2273, (-3, 12, 2.0006, 17.0006))\n",
"\n",
"dt = datetime.datetime(2018, 7, 8)\n",
"start_time = datetime.datetime.combine(dt.date(), datetime.time(0,0))\n",
"end_time = start_time + datetime.timedelta(days=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Find files"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"files = find_files_and_readers(sensor=sensor,\n",
" start_time=start_time,\n",
" end_time=end_time,\n",
" base_dir=base_dir,\n",
" reader=reader)\n",
"\n",
"len(files[reader])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Create Scene"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'bayes_an', 'cloud_an', 'confidence_an', 'latitude', 'longitude', 'pointing_an', 'satellite_azimuth_angle', 'satellite_zenith_angle', 'solar_azimuth_angle', 'solar_zenith_angle']\n"
]
}
],
"source": [
"scn = Scene(filenames=files, reader=reader)\n",
"bands = scn.available_dataset_names()\n",
"\n",
"print(bands)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load (notice the dtype!)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{DatasetID(name='cloud_an', wavelength=None, resolution=500, polarization=None, calibration=None, level=None, modifiers=()): <xarray.DataArray 'cloud_an' (y: 2400, x: 3000)>\n",
" dask.array<shape=(2400, 3000), dtype=uint16, chunksize=(2400, 3000)>\n",
" Dimensions without coordinates: y, x\n",
" Attributes:\n",
" resolution: 500\n",
" wavelength: None\n",
" coordinates: ('longitude', 'latitude')\n",
" level: None\n",
" units: [-]\n",
" _FillValue: 65535\n",
" platform_name: Sentinel-3A\n",
" modifiers: ()\n",
" polarization: None\n",
" file_type: esa_cloudmask_an\n",
" flag_meanings: visible 1.37_threshold 1.6_small_histogram 1.6_larg...\n",
" flag_masks: [ 1 2 4 8 16 32 64 128 ...\n",
" calibration: None\n",
" sensor: slstr\n",
" name: cloud_an\n",
" start_time: 2018-07-08 09:29:39.702487\n",
" end_time: 2018-07-08 09:32:39.393820\n",
" area: Shape: (2400, 3000)\\nLons: <xarray.DataArray 'longi...\n",
" ancillary_variables: []}"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bands = ['cloud_an']\n",
"scn.load(bands)\n",
"scn.datasets"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Resample (notice the dtype!)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Fill value incompatible with integer data using 65535 instead.\n"
]
},
{
"data": {
"text/plain": [
"{DatasetID(name='cloud_an', wavelength=None, resolution=500, polarization=None, calibration=None, level=None, modifiers=()): <xarray.DataArray 'my_index-2b68d27f8ee0e3aabd06539208cc46f0' (y: 2273, x: 2273)>\n",
" dask.array<shape=(2273, 2273), dtype=float64, chunksize=(2273, 2273)>\n",
" Coordinates:\n",
" * y (y) float64 17.0 17.0 17.0 16.99 16.99 16.99 16.99 16.98 16.98 ...\n",
" * x (x) float64 -2.999 -2.997 -2.994 -2.992 -2.99 -2.988 -2.986 ...\n",
" Attributes:\n",
" resolution: 500\n",
" wavelength: None\n",
" coordinates: ('longitude', 'latitude')\n",
" level: None\n",
" units: [-]\n",
" _FillValue: 65535\n",
" platform_name: Sentinel-3A\n",
" modifiers: ()\n",
" polarization: None\n",
" file_type: esa_cloudmask_an\n",
" flag_meanings: visible 1.37_threshold 1.6_small_histogram 1.6_larg...\n",
" flag_masks: [ 1 2 4 8 16 32 64 128 ...\n",
" calibration: None\n",
" sensor: slstr\n",
" name: cloud_an\n",
" start_time: 2018-07-08 09:29:39.702487\n",
" end_time: 2018-07-08 09:32:39.393820\n",
" area: Area ID: tmp\\nDescription: tmp\\nProjection ID: tmp\\...\n",
" ancillary_variables: []}"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scn_warped = scn.resample(destination=area_slstr, resampler='nearest') # no bilinear avail\n",
"scn_warped.datasets"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Save to Geotiff"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Rasterio 1.0+ required for setting colorinterp\n",
"Rasterio 1.0+ required for setting colorinterp\n",
"Rasterio 1.0+ required for setting colorinterp\n"
]
}
],
"source": [
"filename_nodt = os.path.join(base_dir, 'cloud_an_nodt.tif')\n",
"filename_uint16 = os.path.join(base_dir, 'cloud_an_uint16.tif')\n",
"filename_f32 = os.path.join(base_dir, 'cloud_an_float32.tif')\n",
"\n",
"props = dict(writer='GeoTIFF', enhancement_config=False, compress='LZW', tiled='YES')\n",
"\n",
"scn_warped.save_dataset('cloud_an', filename=filename_nodt, **props)\n",
"scn_warped.save_dataset('cloud_an', filename=filename_uint16, dtype=np.uint16, **props)\n",
"scn_warped.save_dataset('cloud_an', filename=filename_f32, dtype=np.float32, **props)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### \"No datatype\" output"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.DataArray (band: 2, y: 2273, x: 2273)>\n",
"[10333058 values with dtype=uint8]\n",
"Coordinates:\n",
" * band (band) int32 1 2\n",
" * y (y) float64 17.0 17.0 17.0 16.99 16.99 16.99 16.99 16.98 16.98 ...\n",
" * x (x) float64 -2.999 -2.997 -2.994 -2.992 -2.99 -2.988 -2.986 ...\n",
"Attributes:\n",
" is_tiled: 1\n",
" transform: (-3.0, 0.0022, 0.0, 17.0006, 0.0, -0.0021999999999999993)\n",
" res: (0.0022, 0.0021999999999999993)\n",
" nodatavals: (nan, nan)\n",
" crs: +init=epsg:4326"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# notice it got two bands\n",
"ds = xr.open_rasterio(filename_nodt)\n",
"ds"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0, 255], dtype=uint8)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = np.asarray(ds.data)\n",
"data = data[~np.isnan(data)]\n",
"np.unique(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Uint16 output"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.DataArray (band: 2, y: 2273, x: 2273)>\n",
"[10333058 values with dtype=uint16]\n",
"Coordinates:\n",
" * band (band) int32 1 2\n",
" * y (y) float64 17.0 17.0 17.0 16.99 16.99 16.99 16.99 16.98 16.98 ...\n",
" * x (x) float64 -2.999 -2.997 -2.994 -2.992 -2.99 -2.988 -2.986 ...\n",
"Attributes:\n",
" is_tiled: 1\n",
" transform: (-3.0, 0.0022, 0.0, 17.0006, 0.0, -0.0021999999999999993)\n",
" res: (0.0022, 0.0021999999999999993)\n",
" nodatavals: (nan, nan)\n",
" crs: +init=epsg:4326"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# notice it got has bands\n",
"ds = xr.open_rasterio(filename_uint16)\n",
"ds"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0, 65535], dtype=uint16)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = np.asarray(ds.data)\n",
"data = data[~np.isnan(data)]\n",
"np.unique(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Float32 output"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.DataArray (band: 1, y: 2273, x: 2273)>\n",
"[5166529 values with dtype=float32]\n",
"Coordinates:\n",
" * band (band) int32 1\n",
" * y (y) float64 17.0 17.0 17.0 16.99 16.99 16.99 16.99 16.98 16.98 ...\n",
" * x (x) float64 -2.999 -2.997 -2.994 -2.992 -2.99 -2.988 -2.986 ...\n",
"Attributes:\n",
" is_tiled: 1\n",
" transform: (-3.0, 0.0022, 0.0, 17.0006, 0.0, -0.0021999999999999993)\n",
" res: (0.0022, 0.0021999999999999993)\n",
" nodatavals: (nan,)\n",
" crs: +init=epsg:4326"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ds = xr.open_rasterio(filename_f32)\n",
"ds"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Miniconda3\\envs\\satpytest\\lib\\site-packages\\xarray\\backends\\rasterio_.py:108: RuntimeWarning: invalid value encountered in greater\n",
" out = self.riods.value.read(band_key, window=tuple(window))\n",
"C:\\Miniconda3\\envs\\satpytest\\lib\\site-packages\\xarray\\backends\\rasterio_.py:108: RuntimeWarning: invalid value encountered in less\n",
" out = self.riods.value.read(band_key, window=tuple(window))\n"
]
},
{
"data": {
"text/plain": [
"array([ 0., 1., 2., 3., 128., 129., 130., 131., 256., 257., 258.,\n",
" 259., 384., 385., 386., 387.], dtype=float32)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data = np.asarray(ds.data)\n",
"data = data[~np.isnan(data)]\n",
"np.unique(data)"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"# bitflag interpretation:\n",
"\n",
"[(1, 'visible'),\n",
" (2, '1.37_threshold'),\n",
" (4, '1.6_small_histogram'),\n",
" (8, '1.6_large_histogram'),\n",
" (16, '2.25_small_histogram'),\n",
" (32, '2.25_large_histogram'),\n",
" (64, '11_spatial_coherence'),\n",
" (128, 'gross_cloud'),\n",
" (256, 'thin_cirrus'),\n",
" (512, 'medium_high'),\n",
" (1024, 'fog_low_stratus'),\n",
" (2048, '11_12_view_difference'),\n",
" (4096, '3.7_11_view_difference'),\n",
" (8192, 'thermal_histogram'),\n",
" (16384, 'spare'),\n",
" (32768, 'spare')]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.5.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment