Created
July 16, 2018 09:46
-
-
Save RutgerK/3663fea47c03013b11a69512a3971419 to your computer and use it in GitHub Desktop.
SLSTR_cloudmask
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"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