Commit 368858e3 authored by Yannick TANGUY's avatar Yannick TANGUY

commit last modifications from David and internal info in readme

parent c01d59a2
......@@ -13,7 +13,11 @@
- Install ipyleaflet : pip install --user ipyleaflet +
- OTB : ./OTB-6.6.1-Linux64.run --target /home/otb/OTB
- Jupyter :
6. Add dataset
6. Configure network settings (proxy & certificates) for CNES environment.
Instructions from : https://gitlab.cnes.fr/hpc/wikiHPC/wikis/connexion-proxy
and https://gitlab.cnes.fr/hpc/wikiHPC/wikis/Acc%C3%A8s-%C3%A0-un-site-https-hors-master-CNES
7. Add dataset
**Note:** Python3 workaround (https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb/issues/1540)
......
......@@ -30,9 +30,7 @@
"from rasterio.plot import show\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# LUT: color_name (string) <-> band (int)\n",
"COLORS_BANDS = {'Red':4, 'Green':3, 'Blue':2, 'Nir':8}\n",
"import os\n",
"\n",
"# Data directory\n",
"DATA_DIR = \"data\"\n",
......@@ -42,9 +40,11 @@
"LEVEL = \"L2A\"\n",
"DATE = \"20180805\"\n",
"\n",
"# Get Band filename\n",
"S2DIR = glob(os.path.join(DATA_DIR, \"SENTINEL2?_{0}-*_{1}_T{2}_*\".format(DATE, LEVEL, MGRS)))[0]\n",
"BAND_FILENAME = os.path.join(S2DIR, os.path.basename(S2DIR) + \"_FRE_B%d.tif\")\n",
"# Color bands\n",
"COLORS_BANDS = ['Blue', 'Green', 'Red', 'NIR']\n",
"\n",
"# Get S2 product\n",
"S2_product = glob(os.path.join(DATA_DIR, \"SENTINEL2?_{0}-*_{1}_T{2}_*\".format(DATE, LEVEL, MGRS)))[0]\n",
"\n",
"# Normalize bands into 0.0 - 1.0 scale\n",
"def normalize(array, quantile=2):\n",
......@@ -54,19 +54,19 @@
" normalized[normalized<0] = 0\n",
" return normalized\n",
"\n",
"# Get visible bands from SENTINEL-2 product\n",
"s2_prod, rasters, quicklooks = {}, {}, {}\n",
"# Plot bands from SENTINEL-2 product\n",
"raster = rasterio.open(S2_product)\n",
"quicklooks = {}\n",
"\n",
"fig, ax = plt.subplots(1,len(COLORS_BANDS)+1, figsize=(21,7))\n",
"for nfig, color in enumerate(COLORS_BANDS):\n",
" s2_prod[color] = BAND_FILENAME % COLORS_BANDS[color]\n",
" rasters[color] = rasterio.open(s2_prod[color])\n",
" quicklooks[color] = normalize(rasters[color].read(1, out_shape=(int(rasters[color].height / 8), int(rasters[color].width / 8))))\n",
" quicklooks[color] = normalize(raster.read(nfig+1, out_shape=(int(raster.height / 8), int(raster.width / 8))))\n",
" try:\n",
" ax[nfig].imshow(quicklooks[color], cmap=color+'s')\n",
" except ValueError:\n",
" ax[nfig].imshow(quicklooks[color], cmap='Oranges')\n",
" ax[nfig].set_title(color+' channel')\n",
" \n",
"\n",
"ax[len(COLORS_BANDS)].imshow(np.dstack((quicklooks['Red'], quicklooks['Green'], quicklooks['Blue'])))\n",
"ax[len(COLORS_BANDS)].set_title('RGB')\n",
"\n",
......@@ -82,20 +82,27 @@
"outputs": [],
"source": [
"from rasterio.warp import transform_bounds\n",
"from ipyleaflet import Map, Rectangle, ImageOverlay, FullScreenControl, DrawControl\n",
"from ipyleaflet import Map, Rectangle, ImageOverlay, FullScreenControl, DrawControl, LayersControl\n",
"\n",
"# Plot footprint on a map\n",
"\n",
"# - get bounding box\n",
"epsg4326 = {'init': 'EPSG:4326'}\n",
"bounds = transform_bounds(rasters['Red'].crs, epsg4326, *rasters['Red'].bounds)\n",
"bounds = transform_bounds(raster.crs, epsg4326, *raster.bounds)\n",
"center = [(bounds[0]+bounds[2])/2, (bounds[1]+bounds[3])/2]\n",
"\n",
"# - get centered map\n",
"m = Map(center=(center[-1], center[0]), zoom=7)\n",
"m = Map(center=(center[-1], center[0]), zoom=10)\n",
"\n",
"# - create quicklook\n",
"import otbApplication\n",
"quicklook_url = os.path.splitext(S2_product)[0] + \"_PREVIEW.JPG\"\n",
"app = otbApplication.Registry.CreateApplication(\"DynamicConvert\")\n",
"app.SetParameterString('in', S2_product)\n",
"app.SetParameterString('out', quicklook_url)\n",
"app.ExecuteAndWriteOutput()\n",
"\n",
"# - plot quicklook\n",
"quicklook_url = os.path.join(S2DIR, os.path.basename(S2DIR) + \"_QKL_ALL.jpg\")\n",
"quicklook = ImageOverlay(\n",
" url=quicklook_url,\n",
" bounds=((bounds[1], bounds[0]),(bounds[3], bounds[2])),\n",
......@@ -127,12 +134,12 @@
" lon_min_max = [np.amin(coordinates[:,0]), np.amax(coordinates[:,0])]\n",
" lat_min_max = [np.amin(coordinates[:,1]), np.amax(coordinates[:,1])]\n",
" lons, lats = np.meshgrid(lon_min_max, lat_min_max)\n",
" xs, ys = warp.transform({'init': 'EPSG:4326'}, rasters['Red'].crs, lons.flatten(), lats.flatten())\n",
" xs, ys = transform({'init': 'EPSG:4326'}, raster.crs, lons.flatten(), lats.flatten())\n",
" \n",
" # Get the region of interest in the image\n",
" rows, cols = [], []\n",
" for x, y in zip(xs, ys):\n",
" row, col = ~rasters['Red'].affine * (x, y)\n",
" row, col = ~raster.affine * (x, y)\n",
" rows.append(row)\n",
" cols.append(col)\n",
" row_min, row_max = np.amin(rows), np.amax(rows)\n",
......@@ -142,11 +149,7 @@
" return startx, starty, endx, endy\n",
"\n",
"startx, starty, endx, endy = get_bounding_box_from_draw()\n",
"red_ext = normalize(rasters['Red'].read(1, window=((startx, endx), (starty, endy))))\n",
"green_ext = normalize(rasters['Green'].read(1, window=((startx, endx), (starty, endy))))\n",
"blue_ext = normalize(rasters['Blue'].read(1, window=((startx, endx), (starty, endy))))\n",
"\n",
"plt.imshow(np.dstack((red_ext, green_ext, blue_ext)))"
"print (startx, starty, endx, endy)"
]
},
{
......@@ -159,9 +162,9 @@
],
"metadata": {
"kernelspec": {
"display_name": "tensorflow",
"display_name": "Python 3",
"language": "python",
"name": "tf"
"name": "python3"
},
"language_info": {
"codemirror_mode": {
......@@ -173,7 +176,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
"version": "3.5.6"
}
},
"nbformat": 4,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment