{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "loading = True" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "3776c009900049449f4ec7be408d3ede", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Accordion(children=(VBox(children=(Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', '…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from utils.utils import (\n", " lat_lon_to_epsg)\n", "import uuid\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "from datacube import Datacube\n", "import xarray\n", "import ipywidgets as widgets\n", "from ipywidgets import Button, Layout,HBox, VBox,Accordion,ToggleButtons,SelectionRangeSlider,Label\n", "import json\n", "import os\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "from traitlets import link\n", "import datetime\n", "from datacube.storage import masking # Import masking capabilities\n", "import pandas as pd\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "dc = Datacube()\n", "\n", "###need to add lat long and address search\n", "\n", "if 'center' not in locals():\n", " center =[-34.42989,116.63979]\n", " \n", "zoom = 2\n", "\n", "product_summary = widgets.Output(layout={'border': '1px solid black'})\n", "#dimensions = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "\n", "def generate_available_datasets_table():\n", " products_from_cube = dc.list_products()\n", " product_list = list(dc.list_products()['name'].values)\n", " product_description = list(dc.list_products()['description'].values)\n", " \n", " build_table = []\n", " \n", " config = load_config('./configIndex.txt')\n", " \n", " for num, product in enumerate(product_list):\n", " ds = dc.find_datasets(product=product,\n", " latitude=tuple(config['lat']),\n", " longitude=tuple(config['lon']),\n", " time=tuple(sorted(config['time']))\n", " )\n", " captured_dates = []\n", " for dataset in ds:\n", " captured_dates.append(dataset.time.begin.date())\n", " if len(captured_dates)>0:\n", " number_epochs = len(set(captured_dates))\n", " number_tiles = len(captured_dates)\n", "\n", " captured_dates_sorted = sorted(captured_dates)\n", "\n", " start_date = captured_dates_sorted[0]\n", " end_date = captured_dates_sorted[len(captured_dates_sorted)-1]\n", "\n", "\n", " build_table.append([product,\n", " product_description[num],\n", " number_tiles,\n", " number_epochs,\n", " start_date,\n", " end_date\n", " ])\n", " build_table = pd.DataFrame(build_table,columns = ['Product Name', \n", " 'Product Description',\n", " 'Number of Tiles',\n", " 'Number of Epochs',\n", " 'Start Date',\n", " 'End Date'])\n", " product_summary.clear_output()\n", " with product_summary:\n", " display(build_table)\n", "\n", "\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", " \n", "def update_config(filename,variable,value):\n", " \n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " data[variable] = value\n", " \n", " os.remove(filename)\n", " with open(filename, 'w') as f:\n", " json.dump(data, f, indent=4)\n", "\n", "m = Map(center=center, zoom=zoom)\n", "m2 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery,layout=m.layout)\n", "\n", "draw_control = DrawControl(rectangle={'shapeOptions': {'color': '#0000FF'}})\n", "draw_control2 = DrawControl(rectangle={'shapeOptions': {'color': '#0000FF'}})\n", "\n", "def handle_draw(self, action, geo_json):\n", " if action == 'created':\n", " m2.add_layer(GeoJSON(data=draw_control.last_draw))\n", " draw_control2.last_draw =draw_control.last_draw\n", " \n", " lon_max = max([draw_control.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][0]])\n", " lon_min = min([draw_control.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][0]])\n", "\n", " lat_max = max([draw_control.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][1]])\n", " lat_min = min([draw_control.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][1]])\n", " EPSG = lat_lon_to_epsg(lat_max,lon_min)\n", " \n", " \n", " \n", " #lat = {'lat' : (lat_max,lat_min),\n", " # 'lon' : (lon_max,lon_min)}\n", " update_config('./configIndex.txt',\n", " 'output_crs',\n", " 'epsg:' + EPSG)\n", " \n", " update_config('./configIndex.txt',\n", " 'lat',\n", " (lat_max,lat_min))\n", "\n", " update_config('./configIndex.txt',\n", " 'lon',\n", " (lon_max,lon_min))\n", " update_config('./configIndex.txt',\n", " 'geoJSON',\n", " draw_control.last_draw\n", " )\n", " update_config('./configIndex.txt',\n", " 'load_id',\n", " str(uuid.uuid4())\n", " )\n", "\n", "\n", " if action == 'deleted':\n", " while len(m2.layers)>1:\n", " m2.remove_layer(m2.layers[1])\n", "\n", "def handle_draw2(self, action, geo_json):\n", " \n", " if action == 'created':\n", " m.add_layer(GeoJSON(data=draw_control2.last_draw))\n", " draw_control.last_draw =draw_control2.last_draw\n", " \n", " lon_max = max([draw_control2.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][0]])\n", " lon_min = min([draw_control2.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][0]])\n", "\n", " lat_max = max([draw_control2.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][1]])\n", " lat_min = min([draw_control2.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][1]])\n", " EPSG = lat_lon_to_epsg(lat_max,lon_min)\n", " \n", " update_config('./configIndex.txt',\n", " 'output_crs',\n", " 'epsg:' + EPSG)\n", " \n", " update_config('./configIndex.txt',\n", " 'lat',\n", " (lat_max,lat_min))\n", "\n", " update_config('./configIndex.txt',\n", " 'lon',\n", " (lon_max,lon_min))\n", " \n", " update_config('./configIndex.txt',\n", " 'geoJSON',\n", " draw_control2.last_draw\n", " )\n", " update_config('./configIndex.txt',\n", " 'load_id',\n", " str(uuid.uuid4())\n", " )\n", "\n", "\n", " if action == 'deleted':\n", " while len(m.layers)>1:\n", " m.remove_layer(m.layers[1])\n", " \n", "#add handlers to draw controls \n", "draw_control.on_draw(handle_draw)\n", "draw_control2.on_draw(handle_draw2)\n", "\n", "#add draw controls to maps\n", "m.add_control(draw_control)\n", "m2.add_control(draw_control2)\n", "\n", "#We can use link to synchronize traitlets of the two maps:\n", "\n", "map_center_link = link((m, 'center'), (m2, 'center'))\n", "map_zoom_link = link((m, 'zoom'), (m2, 'zoom'))\n", "\n", "dates = [datetime.date(1986,1,1) + datetime.timedelta(days =i) for i in range(1,12550)]\n", "\n", "date_range = SelectionRangeSlider(options=dates,\n", " index=(1,12500),\n", " description = 'Date Range',\n", " disabled=False,\n", " layout = Layout(width='100%',height = '100px'))\n", "\n", "def date_func(b):\n", " start_date,end_date = date_range.value\n", " update_config('./configIndex.txt',\n", " 'time',\n", " (start_date.strftime(\"%Y-%m-%d\"),\n", " end_date.strftime(\"%Y-%m-%d\")))\n", "\n", " generate_available_datasets_table()\n", "\n", "button = widgets.Button(description=\"Query Cube\",\n", " button_Style = 'success',\n", " layout=Layout(width='100%', height = '40px'))\n", "\n", "button.on_click(date_func)\n", " \n", "case_study_select = VBox([m,m2,product_summary, date_range, button])\n", " \n", "#date_range.observe(date_func, names = 'value')\n", "\n", "exe_load = True\n", "\n", "accordion = Accordion(children=[case_study_select])\n", "accordion.set_title(0, 'Select Case Study')\n", "accordion" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Macro `select_case_study_area` created. To execute, type its name (without quotes).\n", "=== Macro contents: ===\n", "from utils.utils import (\n", " lat_lon_to_epsg)\n", "import uuid\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "from datacube import Datacube\n", "import xarray\n", "import ipywidgets as widgets\n", "from ipywidgets import Button, Layout,HBox, VBox,Accordion,ToggleButtons,SelectionRangeSlider,Label\n", "import json\n", "import os\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "from traitlets import link\n", "import datetime\n", "from datacube.storage import masking # Import masking capabilities\n", "import pandas as pd\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "dc = Datacube()\n", "\n", "###need to add lat long and address search\n", "\n", "if 'center' not in locals():\n", " center =[-34.42989,116.63979]\n", " \n", "zoom = 2\n", "\n", "product_summary = widgets.Output(layout={'border': '1px solid black'})\n", "#dimensions = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "\n", "def generate_available_datasets_table():\n", " products_from_cube = dc.list_products()\n", " product_list = list(dc.list_products()['name'].values)\n", " product_description = list(dc.list_products()['description'].values)\n", " \n", " build_table = []\n", " \n", " config = load_config('./configIndex.txt')\n", " \n", " for num, product in enumerate(product_list):\n", " ds = dc.find_datasets(product=product,\n", " latitude=tuple(config['lat']),\n", " longitude=tuple(config['lon']),\n", " time=tuple(sorted(config['time']))\n", " )\n", " captured_dates = []\n", " for dataset in ds:\n", " captured_dates.append(dataset.time.begin.date())\n", " if len(captured_dates)>0:\n", " number_epochs = len(set(captured_dates))\n", " number_tiles = len(captured_dates)\n", "\n", " captured_dates_sorted = sorted(captured_dates)\n", "\n", " start_date = captured_dates_sorted[0]\n", " end_date = captured_dates_sorted[len(captured_dates_sorted)-1]\n", "\n", "\n", " build_table.append([product,\n", " product_description[num],\n", " number_tiles,\n", " number_epochs,\n", " start_date,\n", " end_date\n", " ])\n", " build_table = pd.DataFrame(build_table,columns = ['Product Name', \n", " 'Product Description',\n", " 'Number of Tiles',\n", " 'Number of Epochs',\n", " 'Start Date',\n", " 'End Date'])\n", " product_summary.clear_output()\n", " with product_summary:\n", " display(build_table)\n", "\n", "\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", " \n", "def update_config(filename,variable,value):\n", " \n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " data[variable] = value\n", " \n", " os.remove(filename)\n", " with open(filename, 'w') as f:\n", " json.dump(data, f, indent=4)\n", "\n", "m = Map(center=center, zoom=zoom)\n", "m2 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery,layout=m.layout)\n", "\n", "draw_control = DrawControl(rectangle={'shapeOptions': {'color': '#0000FF'}})\n", "draw_control2 = DrawControl(rectangle={'shapeOptions': {'color': '#0000FF'}})\n", "\n", "def handle_draw(self, action, geo_json):\n", " if action == 'created':\n", " m2.add_layer(GeoJSON(data=draw_control.last_draw))\n", " draw_control2.last_draw =draw_control.last_draw\n", " \n", " lon_max = max([draw_control.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][0]])\n", " lon_min = min([draw_control.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][0]])\n", "\n", " lat_max = max([draw_control.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][1]])\n", " lat_min = min([draw_control.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control.last_draw['geometry']['coordinates'][0][2][1]])\n", " EPSG = lat_lon_to_epsg(lat_max,lon_min)\n", " \n", " \n", " \n", " #lat = {'lat' : (lat_max,lat_min),\n", " # 'lon' : (lon_max,lon_min)}\n", " update_config('./configIndex.txt',\n", " 'output_crs',\n", " 'epsg:' + EPSG)\n", " \n", " update_config('./configIndex.txt',\n", " 'lat',\n", " (lat_max,lat_min))\n", "\n", " update_config('./configIndex.txt',\n", " 'lon',\n", " (lon_max,lon_min))\n", " update_config('./configIndex.txt',\n", " 'geoJSON',\n", " draw_control.last_draw\n", " )\n", " update_config('./configIndex.txt',\n", " 'load_id',\n", " str(uuid.uuid4())\n", " )\n", "\n", "\n", " if action == 'deleted':\n", " while len(m2.layers)>1:\n", " m2.remove_layer(m2.layers[1])\n", "\n", "def handle_draw2(self, action, geo_json):\n", " \n", " if action == 'created':\n", " m.add_layer(GeoJSON(data=draw_control2.last_draw))\n", " draw_control.last_draw =draw_control2.last_draw\n", " \n", " lon_max = max([draw_control2.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][0]])\n", " lon_min = min([draw_control2.last_draw['geometry']['coordinates'][0][0][0],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][0]])\n", "\n", " lat_max = max([draw_control2.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][1]])\n", " lat_min = min([draw_control2.last_draw['geometry']['coordinates'][0][0][1],\n", " draw_control2.last_draw['geometry']['coordinates'][0][2][1]])\n", " EPSG = lat_lon_to_epsg(lat_max,lon_min)\n", " \n", " update_config('./configIndex.txt',\n", " 'output_crs',\n", " 'epsg:' + EPSG)\n", " \n", " update_config('./configIndex.txt',\n", " 'lat',\n", " (lat_max,lat_min))\n", "\n", " update_config('./configIndex.txt',\n", " 'lon',\n", " (lon_max,lon_min))\n", " \n", " update_config('./configIndex.txt',\n", " 'geoJSON',\n", " draw_control2.last_draw\n", " )\n", " update_config('./configIndex.txt',\n", " 'load_id',\n", " str(uuid.uuid4())\n", " )\n", "\n", "\n", " if action == 'deleted':\n", " while len(m.layers)>1:\n", " m.remove_layer(m.layers[1])\n", " \n", "#add handlers to draw controls \n", "draw_control.on_draw(handle_draw)\n", "draw_control2.on_draw(handle_draw2)\n", "\n", "#add draw controls to maps\n", "m.add_control(draw_control)\n", "m2.add_control(draw_control2)\n", "\n", "#We can use link to synchronize traitlets of the two maps:\n", "\n", "map_center_link = link((m, 'center'), (m2, 'center'))\n", "map_zoom_link = link((m, 'zoom'), (m2, 'zoom'))\n", "\n", "dates = [datetime.date(1986,1,1) + datetime.timedelta(days =i) for i in range(1,12550)]\n", "\n", "date_range = SelectionRangeSlider(options=dates,\n", " index=(1,12500),\n", " description = 'Date Range',\n", " disabled=False,\n", " layout = Layout(width='100%',height = '100px'))\n", "\n", "def date_func(b):\n", " start_date,end_date = date_range.value\n", " update_config('./configIndex.txt',\n", " 'time',\n", " (start_date.strftime(\"%Y-%m-%d\"),\n", " end_date.strftime(\"%Y-%m-%d\")))\n", "\n", " generate_available_datasets_table()\n", "\n", "button = widgets.Button(description=\"Query Cube\",\n", " button_Style = 'success',\n", " layout=Layout(width='100%', height = '40px'))\n", "\n", "button.on_click(date_func)\n", " \n", "case_study_select = VBox([m,m2,product_summary, date_range, button])\n", " \n", "#date_range.observe(date_func, names = 'value')\n", "\n", "exe_load = True\n", "\n", "accordion = Accordion(children=[case_study_select])\n", "accordion.set_title(0, 'Select Case Study')\n", "accordion\n", " " ] } ], "source": [ "%macro select_case_study_area 2" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stored 'select_case_study_area' (Macro)\n" ] } ], "source": [ "%store select_case_study_area" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n", "No Case Study Selected\n", "No index selected\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "cf18410630154823b1e62dbe51e7967c", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Output(layout=Layout(border='1px solid black')), ToggleButtons(description='Product:', options=…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "##Landsat 8 AWS\n", "%pylab notebook\n", "import sys\n", "import os\n", "import osr \n", "import ogr\n", "import ipywidgets as widgets\n", "import pandas as pd\n", "import rasterio.features\n", "import ipywidgets as widgets\n", "import xarray\n", "#sys.path.append(os.path.expanduser('dea-notebooks/Scripts'))\n", "from utils import BandIndices\n", "from utils.utils import (transform_from_wgs,\n", " transform_from_wgs_poly)\n", "import json\n", "from datacube import Datacube\n", "from datacube.storage import masking # Import masking capabilities\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "if 'case_study' not in locals():\n", " global case_study\n", " case_study = './configIndex.txt'\n", "\n", "\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "\n", "global output_pandas\n", "output_pandas= pd.DataFrame() \n", "\n", "global geoJSONs\n", "geoJSONs = []\n", "\n", "global datasetID\n", "datasetID = 0\n", "\n", "###Load EPSG\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", "\n", "##Define output widget to display table\n", "out = widgets.Output(layout={'border': '1px solid black'})\n", "info = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "def load_into_cube_aws_landsat8(measurements = None, cloud_mask = False):\n", "\n", " def lat_lon_to_epsg(lat_max, lon_min):\n", " return str(int(32700 - round((45 + lat_max) / 90, 0) * 100 + round((183 + lon_min) / 6, 0)))\n", "\n", " \n", " data = load_config(case_study)\n", "\n", " \n", " lats= data['lat']\n", " longs= data['lon']\n", " global EPSG\n", " EPSG = lat_lon_to_epsg(lats[1],longs[0])\n", " \n", " ds = dc.load(product='ls8_level1_usgs', \n", " measurements = ('red','green','blue','nir','swir1','swir2','quality'),\n", " group_by='solar_day',\n", " y=tuple(data['lat']),\n", " x=tuple(data['lon']),\n", " time=tuple(data['time']),\n", " output_crs='epsg:' + EPSG,\n", " resolution = (-25,25))\n", "\n", " #filter for clouds\n", " clean_pixel_mask = masking.make_mask(\n", " ds.quality,\n", " cloud=False,\n", " radiometric_saturation='none',\n", " terrain_occlusion = False)\n", "\n", " masked_cloud = ds.where(clean_pixel_mask)\n", "\n", " return masked_cloud\n", "\n", "#### Get geometry from config file draw polygon extent.\n", "config =load_config(case_study)\n", "\n", "geoJSON_Extent = config['geoJSON']\n", "#m3.add_layer(GeoJSON(data=draw_control2.last_draw))\n", "\n", "info_last_geom = ogr.CreateGeometryFromJson(str(geoJSON_Extent['geometry']))\n", "if info_last_geom:\n", " zoom = 13\n", " center = [info_last_geom.Centroid().GetY(),info_last_geom.Centroid().GetX()]\n", " m3 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "else:\n", " m3 = Map(center=center, zoom=zoom)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "\n", "dc = Datacube()\n", "\n", "if 'loading' in locals():\n", " ds = None\n", " print ('No Case Study Selected')\n", " target_dataset,desc = None, 'No index selected'\n", "else:\n", " geoJSON_Extent['properties']['style']['fillOpacity'] = 0\n", " geoJSON_Extent['properties']['style']['color'] = 'red'\n", " \n", " if 'loaded' in locals():\n", " if loaded != config['load_id']:\n", " print ('Retrieving New Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " else:\n", " print ('Retrieving Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " \n", " target_dataset,desc = BandIndices.calculate_indices(ds,'NDVI')\n", "\n", " m3.add_layer(GeoJSON(data=geoJSON_Extent))\n", "\n", "\n", "\n", "indicesa = ['NDVI','GNDVI','NDWI', 'NDMI', 'NDBI', 'NBR']\n", "\n", "with info:\n", " print(desc)\n", "\n", "indices_buttons = widgets.ToggleButtons(\n", " options=indicesa,\n", " description='Product:',\n", " disabled=False,\n", " button_style='' \n", " #tooltips=indices_description,\n", " )\n", "\n", "def indices_func(change):\n", " global target_dataset\n", " target_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", "\n", " info.clear_output()\n", " with info:\n", " print(desc)\n", "\n", "indices_buttons.observe(indices_func, names = 'value')\n", "\n", "\n", "dataP = load_config(case_study)\n", "\n", "#EPSG = int(dataP['output_crs'].split(':')[1])\n", "#EPSG = 3577\n", "labels = []\n", "def plot_nDVI(self, action, geo_json):\n", " global output_pandas\n", " global datasetID\n", "\n", " plt.figure(0,[10,5])\n", " plt.ylim(-1, 1)\n", "\n", " #print (geo_json)\n", " if geo_json['geometry']['type'] == 'Point':\n", "\n", " pixel_drill_long,pixel_drill_lat = geo_json['geometry']['coordinates']\n", "\n", " pixel_drill_x, pixel_drill_y = transform_from_wgs(pixel_drill_long,pixel_drill_lat,int(EPSG))\n", " drill_cube = target_dataset.sel(x=[pixel_drill_x,pixel_drill_x + 25], y=[pixel_drill_y, pixel_drill_y + 25],\n", " method = 'nearest')\n", "\n", " time = np.ravel(drill_cube.isel(x=[0],y=[0]).time)\n", " values = np.ravel(drill_cube.isel(x=[0],y=[0]).values)\n", "\n", " #drill_cube.isel(x=[0],y=[0]).plot()\n", " label = str(indices_buttons.value) + ' point ' + str(datasetID)\n", " labels.append(label)\n", " \n", " xarray.plot.plot(drill_cube.isel(x=[0],y=[0]).interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", " datasetID = datasetID + 1\n", " \n", "\n", " if geo_json['geometry']['type'] == 'Polygon':\n", "\n", " geom = transform_from_wgs_poly(geo_json['geometry'],int(EPSG))\n", "\n", " mask = rasterio.features.geometry_mask([geom for geoms in [geom]],\n", " out_shape=ds.geobox.shape,\n", " transform=ds.geobox.affine,\n", " all_touched=False,\n", " invert=True)\n", " dataMasked = target_dataset.where(mask)\n", " data_masked_mean = dataMasked.mean(dim = ['x','y'])\n", " #data_masked_mean.plot()\n", "\n", " time = np.ravel(data_masked_mean.time)\n", " values = np.ravel(data_masked_mean.values)\n", "\n", " label = str(indices_buttons.value) + ' polygon ' + str(datasetID)\n", " labels.append(label)\n", "\n", " xarray.plot.plot(data_masked_mean.interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", "\n", " datasetID = datasetID + 1\n", "\n", " geoJSONs.append(geo_json)\n", "\n", " output_pandas = pd.concat([output_pandas,pd.DataFrame(values,index = time, columns = [label])],axis = 1)\n", " output_pandas.dropna(thresh=datasetID,inplace=True)\n", " out.clear_output()\n", " with out:\n", " display(output_pandas)\n", "\n", "draw_control3.on_draw(plot_nDVI)\n", "m3.add_control(draw_control3)\n", "\n", "pixel_drill = widgets.VBox([info,indices_buttons,m3,out])\n", "\n", "pixel_drill" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Macro `band_indices_aws_landsat8` created. To execute, type its name (without quotes).\n", "=== Macro contents: ===\n", "##Landsat 8 AWS\n", "get_ipython().run_line_magic('pylab', 'notebook')\n", "import sys\n", "import os\n", "import osr \n", "import ogr\n", "import ipywidgets as widgets\n", "import pandas as pd\n", "import rasterio.features\n", "import ipywidgets as widgets\n", "import xarray\n", "#sys.path.append(os.path.expanduser('dea-notebooks/Scripts'))\n", "from utils import BandIndices\n", "from utils.utils import (transform_from_wgs,\n", " transform_from_wgs_poly)\n", "import json\n", "from datacube import Datacube\n", "from datacube.storage import masking # Import masking capabilities\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "if 'case_study' not in locals():\n", " global case_study\n", " case_study = './configIndex.txt'\n", "\n", "\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "\n", "global output_pandas\n", "output_pandas= pd.DataFrame() \n", "\n", "global geoJSONs\n", "geoJSONs = []\n", "\n", "global datasetID\n", "datasetID = 0\n", "\n", "###Load EPSG\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", "\n", "##Define output widget to display table\n", "out = widgets.Output(layout={'border': '1px solid black'})\n", "info = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "def load_into_cube_aws_landsat8(measurements = None, cloud_mask = False):\n", "\n", " def lat_lon_to_epsg(lat_max, lon_min):\n", " return str(int(32700 - round((45 + lat_max) / 90, 0) * 100 + round((183 + lon_min) / 6, 0)))\n", "\n", " \n", " data = load_config(case_study)\n", "\n", " \n", " lats= data['lat']\n", " longs= data['lon']\n", " global EPSG\n", " EPSG = lat_lon_to_epsg(lats[1],longs[0])\n", " \n", " ds = dc.load(product='ls8_level1_usgs', \n", " measurements = ('red','green','blue','nir','swir1','swir2','quality'),\n", " group_by='solar_day',\n", " y=tuple(data['lat']),\n", " x=tuple(data['lon']),\n", " time=tuple(data['time']),\n", " output_crs='epsg:' + EPSG,\n", " resolution = (-25,25))\n", "\n", " #filter for clouds\n", " clean_pixel_mask = masking.make_mask(\n", " ds.quality,\n", " cloud=False,\n", " radiometric_saturation='none',\n", " terrain_occlusion = False)\n", "\n", " masked_cloud = ds.where(clean_pixel_mask)\n", "\n", " return masked_cloud\n", "\n", "#### Get geometry from config file draw polygon extent.\n", "config =load_config(case_study)\n", "\n", "geoJSON_Extent = config['geoJSON']\n", "#m3.add_layer(GeoJSON(data=draw_control2.last_draw))\n", "\n", "info_last_geom = ogr.CreateGeometryFromJson(str(geoJSON_Extent['geometry']))\n", "if info_last_geom:\n", " zoom = 13\n", " center = [info_last_geom.Centroid().GetY(),info_last_geom.Centroid().GetX()]\n", " m3 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "else:\n", " m3 = Map(center=center, zoom=zoom)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "\n", "dc = Datacube()\n", "\n", "if 'loading' in locals():\n", " ds = None\n", " print ('No Case Study Selected')\n", " target_dataset,desc = None, 'No index selected'\n", "else:\n", " geoJSON_Extent['properties']['style']['fillOpacity'] = 0\n", " geoJSON_Extent['properties']['style']['color'] = 'red'\n", " \n", " if 'loaded' in locals():\n", " if loaded != config['load_id']:\n", " print ('Retrieving New Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " else:\n", " print ('Retrieving Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " \n", " target_dataset,desc = BandIndices.calculate_indices(ds,'NDVI')\n", "\n", " m3.add_layer(GeoJSON(data=geoJSON_Extent))\n", "\n", "\n", "\n", "indicesa = ['NDVI','GNDVI','NDWI', 'NDMI', 'NDBI', 'NBR']\n", "\n", "with info:\n", " print(desc)\n", "\n", "indices_buttons = widgets.ToggleButtons(\n", " options=indicesa,\n", " description='Product:',\n", " disabled=False,\n", " button_style='' \n", " #tooltips=indices_description,\n", " )\n", "\n", "def indices_func(change):\n", " global target_dataset\n", " target_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", "\n", " info.clear_output()\n", " with info:\n", " print(desc)\n", "\n", "indices_buttons.observe(indices_func, names = 'value')\n", "\n", "\n", "dataP = load_config(case_study)\n", "\n", "#EPSG = int(dataP['output_crs'].split(':')[1])\n", "#EPSG = 3577\n", "labels = []\n", "def plot_nDVI(self, action, geo_json):\n", " global output_pandas\n", " global datasetID\n", "\n", " plt.figure(0,[10,5])\n", " plt.ylim(-1, 1)\n", "\n", " #print (geo_json)\n", " if geo_json['geometry']['type'] == 'Point':\n", "\n", " pixel_drill_long,pixel_drill_lat = geo_json['geometry']['coordinates']\n", "\n", " pixel_drill_x, pixel_drill_y = transform_from_wgs(pixel_drill_long,pixel_drill_lat,int(EPSG))\n", " drill_cube = target_dataset.sel(x=[pixel_drill_x,pixel_drill_x + 25], y=[pixel_drill_y, pixel_drill_y + 25],\n", " method = 'nearest')\n", "\n", " time = np.ravel(drill_cube.isel(x=[0],y=[0]).time)\n", " values = np.ravel(drill_cube.isel(x=[0],y=[0]).values)\n", "\n", " #drill_cube.isel(x=[0],y=[0]).plot()\n", " label = str(indices_buttons.value) + ' point ' + str(datasetID)\n", " labels.append(label)\n", " \n", " xarray.plot.plot(drill_cube.isel(x=[0],y=[0]).interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", " datasetID = datasetID + 1\n", " \n", "\n", " if geo_json['geometry']['type'] == 'Polygon':\n", "\n", " geom = transform_from_wgs_poly(geo_json['geometry'],int(EPSG))\n", "\n", " mask = rasterio.features.geometry_mask([geom for geoms in [geom]],\n", " out_shape=ds.geobox.shape,\n", " transform=ds.geobox.affine,\n", " all_touched=False,\n", " invert=True)\n", " dataMasked = target_dataset.where(mask)\n", " data_masked_mean = dataMasked.mean(dim = ['x','y'])\n", " #data_masked_mean.plot()\n", "\n", " time = np.ravel(data_masked_mean.time)\n", " values = np.ravel(data_masked_mean.values)\n", "\n", " label = str(indices_buttons.value) + ' polygon ' + str(datasetID)\n", " labels.append(label)\n", "\n", " xarray.plot.plot(data_masked_mean.interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", "\n", " datasetID = datasetID + 1\n", "\n", " geoJSONs.append(geo_json)\n", "\n", " output_pandas = pd.concat([output_pandas,pd.DataFrame(values,index = time, columns = [label])],axis = 1)\n", " output_pandas.dropna(thresh=datasetID,inplace=True)\n", " out.clear_output()\n", " with out:\n", " display(output_pandas)\n", "\n", "draw_control3.on_draw(plot_nDVI)\n", "m3.add_control(draw_control3)\n", "\n", "pixel_drill = widgets.VBox([info,indices_buttons,m3,out])\n", "\n", "pixel_drill\n", " " ] } ], "source": [ "%macro band_indices_aws_landsat8 5" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stored 'band_indices_aws_landsat8' (Macro)\n" ] } ], "source": [ "%store band_indices_aws_landsat8" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n", "No Case Study Selected\n", "No index selected\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "51507e92feb9495ca31c33af38399b0d", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Output(layout=Layout(border='1px solid black')), ToggleButtons(description='Product:', options=…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "##Landsat 8 AWS WOFS Time series!!!\n", "%pylab notebook\n", "from importlib import reload # Python 3.4+ only.\n", "import sys\n", "import os\n", "import osr \n", "import ogr\n", "import ipywidgets as widgets\n", "import pandas as pd\n", "import rasterio.features\n", "import ipywidgets as widgets\n", "import xarray\n", "from utils import BandIndices\n", "from matplotlib.colors import from_levels_and_colors\n", "\n", "#sys.path.append(os.path.expanduser('dea-notebooks/Scripts'))\n", "\n", "from utils import Water_Classifier\n", "Water_Classifier = reload(Water_Classifier)\n", "from utils.utils import (transform_from_wgs,\n", " transform_from_wgs_poly)\n", "import json\n", "from datacube import Datacube\n", "from datacube.storage import masking # Import masking capabilities\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "if 'case_study' not in locals():\n", " global case_study\n", " case_study = './configIndex.txt'\n", "\n", "\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "\n", "global output_pandas\n", "output_pandas= pd.DataFrame() \n", "\n", "global geoJSONs\n", "geoJSONs = []\n", "\n", "global datasetID\n", "datasetID = 0\n", "\n", "###Load EPSG\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", "\n", "##Define output widget to display table\n", "out = widgets.Output(layout={'border': '1px solid black'})\n", "info = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "def load_into_cube_aws_landsat8(measurements = None, cloud_mask = False):\n", "\n", " def lat_lon_to_epsg(lat_max, lon_min):\n", " return str(int(32700 - round((45 + lat_max) / 90, 0) * 100 + round((183 + lon_min) / 6, 0)))\n", "\n", " \n", " data = load_config(case_study)\n", "\n", " \n", " lats= data['lat']\n", " longs= data['lon']\n", " global EPSG\n", " EPSG = lat_lon_to_epsg(lats[1],longs[0])\n", " \n", " ds = dc.load(product='ls8_level1_usgs', \n", " measurements = ('red','green','blue','nir','swir1','swir2','quality'),\n", " group_by='solar_day',\n", " y=tuple(data['lat']),\n", " x=tuple(data['lon']),\n", " time=tuple(data['time']),\n", " output_crs='epsg:' + EPSG,\n", " resolution = (-25,25))\n", "\n", " #filter for clouds\n", " clean_pixel_mask = masking.make_mask(\n", " ds.quality,\n", " cloud=False,\n", " radiometric_saturation='none',\n", " terrain_occlusion = False)\n", " masked_cloud = ds.where(clean_pixel_mask)\n", "\n", "\n", " return masked_cloud\n", "\n", "#### Get geometry from config file draw polygon extent.\n", "config =load_config(case_study)\n", "\n", "geoJSON_Extent = config['geoJSON']\n", "#m3.add_layer(GeoJSON(data=draw_control2.last_draw))\n", "\n", "info_last_geom = ogr.CreateGeometryFromJson(str(geoJSON_Extent['geometry']))\n", "if info_last_geom:\n", " zoom = 13\n", " center = [info_last_geom.Centroid().GetY(),info_last_geom.Centroid().GetX()]\n", " m3 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "else:\n", " m3 = Map(center=center, zoom=zoom)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "\n", "dc = Datacube()\n", "\n", "indicesa = ['NDWI', 'WOFS']\n", "\n", "\n", "indices_buttons = widgets.ToggleButtons(\n", " options=indicesa,\n", " description='Product:',\n", " disabled=False,\n", " button_style='' \n", " #tooltips=indices_description,\n", " )\n", "\n", "\n", "if 'loading' in locals():\n", " ds = None\n", " print ('No Case Study Selected')\n", " target_dataset,desc = None, 'No index selected'\n", "else:\n", " geoJSON_Extent['properties']['style']['fillOpacity'] = 0\n", " geoJSON_Extent['properties']['style']['color'] = 'red'\n", " \n", " if 'loaded' in locals():\n", " if loaded != config['load_id']:\n", " print ('Retrieving New Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", " else:\n", " print ('Retrieving Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", " target_dataset = (target_dataset * 625) / 10000\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999)\n", "\n", " # Load data the first time.\n", "# indices_func(True)\n", " \n", " m3.add_layer(GeoJSON(data=geoJSON_Extent))\n", "\n", " \n", "def indices_func(change):\n", " global target_dataset\n", " if indices_buttons.value == 'NDWI':\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", "\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", " target_dataset = (target_dataset * 625) / 10000\n", " else:\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999)\n", " target_dataset,desc = Water_Classifier.water_classifier(dse)\n", " target_dataset = (target_dataset * 625) / 10000\n", "\n", " info.clear_output()\n", " with info:\n", " print(desc)\n", "\n", "\n", "with info:\n", " print(desc)\n", " \n", "indices_buttons.observe(indices_func, names = 'value')\n", "\n", "dataP = load_config(case_study)\n", "\n", "#EPSG = int(dataP['output_crs'].split(':')[1])\n", "#EPSG = 3577\n", "labels = []\n", "\n", "def onclick(event):\n", " global timeOfInterest\n", " timeOfInterest = event.xdata\n", " time_slice = matplotlib.dates.num2date(timeOfInterest).date()\n", " time_slice = str(time_slice)\n", " time_slice = pd.to_datetime(time_slice, format='%Y-%m-%d')\n", " label_html.value = 'Date selected from stack plot is ' + str(time_slice)[:10] + ''\n", "\n", "def display_extent(active):\n", " \n", " time_slice = matplotlib.dates.num2date(timeOfInterest).date()\n", " time_slice = str(time_slice)\n", " time_slice = pd.to_datetime(time_slice, format='%Y-%m-%d')\n", "\n", " plot_graph = cleaned_mask.sel(time=time_slice, method='nearest')\n", " \n", " plot_graph = plot_graph.wofs[firstrow:lastrow,firstcol:lastcol]\n", " cmap, norm = from_levels_and_colors([-10000,0,0.04,2],['red','green','blue'])\n", " \n", " if plt.fignum_exists(collect + 1):\n", " plt.figure(num= collect + 1, figsize=(7.5, 7.5))\n", " plt.cla()\n", "\n", " plt.imshow(plot_graph,cmap=cmap, norm=norm)\n", " else:\n", " plt.figure(num= collect + 1, figsize=(7.5, 7.5))\n", " \n", " plt.imshow(plot_graph,cmap=cmap, norm=norm)\n", " \n", " plt.title( 'Water Extent on ' + str(time_slice)[:10])\n", "\n", " cbar = plt.colorbar(orientation = 'horizontal') \n", "\n", " cbar.ax.get_xaxis().set_ticks([])\n", " for j, lab in enumerate(['no data','land','water']):\n", " cbar.ax.text((2.5 * j + 1) / 8.0, 0.5, lab, ha='center', va='center')\n", " cbar.ax.get_xaxis().labelpad = 10\n", " cbar.ax.set_xlabel('classification')\n", "\n", "\n", " \n", "collect = 0\n", "\n", "def plot_nDVI(self, action, geo_json):\n", " global output_pandas\n", " global datasetID\n", " global mask\n", " global cleaned_mask\n", " global collect\n", " global firstcol,firstrow,lastcol,lastrow\n", " #print (geo_json)\n", " if geo_json['geometry']['type'] == 'Point':\n", "\n", " pixel_drill_long,pixel_drill_lat = geo_json['geometry']['coordinates']\n", "\n", " pixel_drill_x, pixel_drill_y = transform_from_wgs(pixel_drill_long,pixel_drill_lat,int(EPSG))\n", " drill_cube = target_dataset.sel(x=[pixel_drill_x,pixel_drill_x + 25], y=[pixel_drill_y, pixel_drill_y + 25],\n", " method = 'nearest')\n", "\n", " time = np.ravel(drill_cube.isel(x=[0],y=[0]).time)\n", " values = np.ravel(drill_cube.isel(x=[0],y=[0]).values)\n", "\n", " #drill_cube.isel(x=[0],y=[0]).plot()\n", " label = str(indices_buttons.value) + ' point ' + str(datasetID)\n", " labels.append(label)\n", " \n", " xarray.plot.plot(drill_cube.isel(x=[0],y=[0]).interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", " datasetID = datasetID + 1\n", " \n", "\n", " if geo_json['geometry']['type'] == 'Polygon':\n", "\n", " geom = transform_from_wgs_poly(geo_json['geometry'],int(EPSG))\n", "\n", " mask = rasterio.features.geometry_mask([geom for geoms in [geom]],\n", " out_shape=ds.geobox.shape,\n", " transform=ds.geobox.affine,\n", " all_touched=False,\n", " invert=True)\n", " \n", " dataMasked = target_dataset.where(mask)\n", " #data_masked_max = dataMasked.resample(time='1M').max()\n", " \n", " orig = dse.where(mask)\n", " #orig_max = orig.resample(time='1M').median()\n", " \n", " no_data_mask = orig.red.where(orig.red==-9999)\n", " no_data = no_data_mask.count(dim = ['x','y'])\n", " \n", " ## Prepare data for ploting extents\n", " cleaned_mask = no_data_mask.combine_first(dataMasked)\n", " \n", " nan_mask = mask.astype('float')\n", " nan_mask[~mask] = np.nan\n", " \n", " nans = np.isnan(nan_mask)\n", " nancols = np.all(nans, axis=0) \n", " nanrows = np.all(nans, axis=1) \n", " \n", " firstcol = nancols.argmin() \n", " firstrow = nanrows.argmin() \n", "\n", " lastcol = len(nancols) - nancols[::-1].argmin() # 8, last index where not NAN\n", " lastrow = len(nanrows) - nanrows[::-1].argmin() # 10\n", "\n", " # Area in Hectares\n", " no_data = (no_data * 625) / 10000\n", " total_area = (mask.sum() * 625) /10000\n", " \n", " data_masked = dataMasked.sum(dim = ['x','y'])\n", " \n", " other_masked = total_area - (data_masked + no_data)\n", "\n", " label = str(indices_buttons.value) + ' polygon ' + str(datasetID)\n", " labels.append(label)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", "\n", "\n", " geoJSONs.append(geo_json)\n", "\n", " fig = plt.figure(num=collect + 1, figsize=(10, 3.5))\n", " time_from = np.ravel(data_masked.time)\n", " \n", " plt.ylabel('Area (Ha)')\n", " color_scheme = '#1f77b4','#2ca02c', '#ff0000' ,'#d62728', \n", " plt.set_title = label\n", " #plt.set_yscale('log')\n", " plt.stackplot(time_from,\n", " data_masked.wofs.values, \n", " other_masked.wofs.values, \n", " no_data.values, \n", " colors = color_scheme,\n", " labels=['Water','Dry','No Data'])\n", " plt.legend(loc='upper left')\n", " fig.canvas.mpl_connect('button_press_event', onclick)\n", " collect = collect + 1 \n", " \n", " output_pandas = pd.concat([output_pandas,pd.DataFrame(data_masked.wofs.values,index = time_from, columns = [label])],axis = 1)\n", " #output_pandas.dropna(thresh=datasetID,inplace=True)\n", "\n", " out.clear_output()\n", " with out:\n", " ##Calc Stats.\n", " display(\"Mean Amount of Area (Ha) Classified as water over time period\")\n", " display(output_pandas.mean())\n", " display(\"Max Amount of Area (Ha) Classified as water over time period\")\n", " display(output_pandas.max())\n", " display(\"Min Amount of Area (Ha) Classified as water over time period\")\n", " display(output_pandas.min())\n", " \n", "\n", "draw_control3.on_draw(plot_nDVI)\n", "m3.add_control(draw_control3)\n", "\n", "label_html = widgets.HTML(\n", " value=\"Date selected from stack plot is \",\n", " layout=widgets.Layout(flex='1 1 auto', width='auto'))\n", "\n", "#box_layout = Layout(display='flex',\n", "# flex_flow='column',\n", "# align_items='stretch',\n", "# border='solid',\n", "# width='50%')\n", "\n", "button_plot = widgets.Button(description='Plot Spatial Extent', \n", " disabled=False,\n", " button_style='', # 'success', 'info', 'warning', 'danger' or ''\n", " tooltip='Click me', layout=widgets.Layout(flex='1 1 auto', width='auto'))\n", "\n", "button_plot.on_click(display_extent)\n", "\n", "\n", "plot_module = widgets.HBox([label_html, button_plot])\n", "\n", "pixel_drill = widgets.VBox([info,indices_buttons,m3,out,plot_module])\n", "\n", "pixel_drill" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Macro `water_stackplot_aws_landsat8` created. To execute, type its name (without quotes).\n", "=== Macro contents: ===\n", "##Landsat 8 AWS WOFS Time series!!!\n", "get_ipython().run_line_magic('pylab', 'notebook')\n", "from importlib import reload # Python 3.4+ only.\n", "import sys\n", "import os\n", "import osr \n", "import ogr\n", "import ipywidgets as widgets\n", "import pandas as pd\n", "import rasterio.features\n", "import ipywidgets as widgets\n", "import xarray\n", "from utils import BandIndices\n", "from matplotlib.colors import from_levels_and_colors\n", "\n", "#sys.path.append(os.path.expanduser('dea-notebooks/Scripts'))\n", "\n", "from utils import Water_Classifier\n", "Water_Classifier = reload(Water_Classifier)\n", "from utils.utils import (transform_from_wgs,\n", " transform_from_wgs_poly)\n", "import json\n", "from datacube import Datacube\n", "from datacube.storage import masking # Import masking capabilities\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "if 'case_study' not in locals():\n", " global case_study\n", " case_study = './configIndex.txt'\n", "\n", "\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "\n", "global output_pandas\n", "output_pandas= pd.DataFrame() \n", "\n", "global geoJSONs\n", "geoJSONs = []\n", "\n", "global datasetID\n", "datasetID = 0\n", "\n", "###Load EPSG\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", "\n", "##Define output widget to display table\n", "out = widgets.Output(layout={'border': '1px solid black'})\n", "info = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "def load_into_cube_aws_landsat8(measurements = None, cloud_mask = False):\n", "\n", " def lat_lon_to_epsg(lat_max, lon_min):\n", " return str(int(32700 - round((45 + lat_max) / 90, 0) * 100 + round((183 + lon_min) / 6, 0)))\n", "\n", " \n", " data = load_config(case_study)\n", "\n", " \n", " lats= data['lat']\n", " longs= data['lon']\n", " global EPSG\n", " EPSG = lat_lon_to_epsg(lats[1],longs[0])\n", " \n", " ds = dc.load(product='ls8_level1_usgs', \n", " measurements = ('red','green','blue','nir','swir1','swir2','quality'),\n", " group_by='solar_day',\n", " y=tuple(data['lat']),\n", " x=tuple(data['lon']),\n", " time=tuple(data['time']),\n", " output_crs='epsg:' + EPSG,\n", " resolution = (-25,25))\n", "\n", " #filter for clouds\n", " clean_pixel_mask = masking.make_mask(\n", " ds.quality,\n", " cloud=False,\n", " radiometric_saturation='none',\n", " terrain_occlusion = False)\n", " masked_cloud = ds.where(clean_pixel_mask)\n", "\n", "\n", " return masked_cloud\n", "\n", "#### Get geometry from config file draw polygon extent.\n", "config =load_config(case_study)\n", "\n", "geoJSON_Extent = config['geoJSON']\n", "#m3.add_layer(GeoJSON(data=draw_control2.last_draw))\n", "\n", "info_last_geom = ogr.CreateGeometryFromJson(str(geoJSON_Extent['geometry']))\n", "if info_last_geom:\n", " zoom = 13\n", " center = [info_last_geom.Centroid().GetY(),info_last_geom.Centroid().GetX()]\n", " m3 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "else:\n", " m3 = Map(center=center, zoom=zoom)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "\n", "dc = Datacube()\n", "\n", "indicesa = ['NDWI', 'WOFS']\n", "\n", "\n", "indices_buttons = widgets.ToggleButtons(\n", " options=indicesa,\n", " description='Product:',\n", " disabled=False,\n", " button_style='' \n", " #tooltips=indices_description,\n", " )\n", "\n", "\n", "if 'loading' in locals():\n", " ds = None\n", " print ('No Case Study Selected')\n", " target_dataset,desc = None, 'No index selected'\n", "else:\n", " geoJSON_Extent['properties']['style']['fillOpacity'] = 0\n", " geoJSON_Extent['properties']['style']['color'] = 'red'\n", " \n", " if 'loaded' in locals():\n", " if loaded != config['load_id']:\n", " print ('Retrieving New Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", " else:\n", " print ('Retrieving Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", " target_dataset = (target_dataset * 625) / 10000\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999)\n", "\n", " # Load data the first time.\n", "# indices_func(True)\n", " \n", " m3.add_layer(GeoJSON(data=geoJSON_Extent))\n", "\n", " \n", "def indices_func(change):\n", " global target_dataset\n", " if indices_buttons.value == 'NDWI':\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", "\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", " target_dataset = (target_dataset * 625) / 10000\n", " else:\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999)\n", " target_dataset,desc = Water_Classifier.water_classifier(dse)\n", " target_dataset = (target_dataset * 625) / 10000\n", "\n", " info.clear_output()\n", " with info:\n", " print(desc)\n", "\n", "\n", "with info:\n", " print(desc)\n", " \n", "indices_buttons.observe(indices_func, names = 'value')\n", "\n", "dataP = load_config(case_study)\n", "\n", "#EPSG = int(dataP['output_crs'].split(':')[1])\n", "#EPSG = 3577\n", "labels = []\n", "\n", "def onclick(event):\n", " global timeOfInterest\n", " timeOfInterest = event.xdata\n", " time_slice = matplotlib.dates.num2date(timeOfInterest).date()\n", " time_slice = str(time_slice)\n", " time_slice = pd.to_datetime(time_slice, format='%Y-%m-%d')\n", " label_html.value = 'Date selected from stack plot is ' + str(time_slice)[:10] + ''\n", "\n", "def display_extent(active):\n", " \n", " time_slice = matplotlib.dates.num2date(timeOfInterest).date()\n", " time_slice = str(time_slice)\n", " time_slice = pd.to_datetime(time_slice, format='%Y-%m-%d')\n", "\n", " plot_graph = cleaned_mask.sel(time=time_slice, method='nearest')\n", " \n", " plot_graph = plot_graph.wofs[firstrow:lastrow,firstcol:lastcol]\n", " cmap, norm = from_levels_and_colors([-10000,0,0.04,2],['red','green','blue'])\n", " \n", " if plt.fignum_exists(collect + 1):\n", " plt.figure(num= collect + 1, figsize=(7.5, 7.5))\n", " plt.cla()\n", "\n", " plt.imshow(plot_graph,cmap=cmap, norm=norm)\n", " else:\n", " plt.figure(num= collect + 1, figsize=(7.5, 7.5))\n", " \n", " plt.imshow(plot_graph,cmap=cmap, norm=norm)\n", " \n", " plt.title( 'Water Extent on ' + str(time_slice)[:10])\n", "\n", " cbar = plt.colorbar(orientation = 'horizontal') \n", "\n", " cbar.ax.get_xaxis().set_ticks([])\n", " for j, lab in enumerate(['no data','land','water']):\n", " cbar.ax.text((2.5 * j + 1) / 8.0, 0.5, lab, ha='center', va='center')\n", " cbar.ax.get_xaxis().labelpad = 10\n", " cbar.ax.set_xlabel('classification')\n", "\n", "\n", " \n", "collect = 0\n", "\n", "def plot_nDVI(self, action, geo_json):\n", " global output_pandas\n", " global datasetID\n", " global mask\n", " global cleaned_mask\n", " global collect\n", " global firstcol,firstrow,lastcol,lastrow\n", " #print (geo_json)\n", " if geo_json['geometry']['type'] == 'Point':\n", "\n", " pixel_drill_long,pixel_drill_lat = geo_json['geometry']['coordinates']\n", "\n", " pixel_drill_x, pixel_drill_y = transform_from_wgs(pixel_drill_long,pixel_drill_lat,int(EPSG))\n", " drill_cube = target_dataset.sel(x=[pixel_drill_x,pixel_drill_x + 25], y=[pixel_drill_y, pixel_drill_y + 25],\n", " method = 'nearest')\n", "\n", " time = np.ravel(drill_cube.isel(x=[0],y=[0]).time)\n", " values = np.ravel(drill_cube.isel(x=[0],y=[0]).values)\n", "\n", " #drill_cube.isel(x=[0],y=[0]).plot()\n", " label = str(indices_buttons.value) + ' point ' + str(datasetID)\n", " labels.append(label)\n", " \n", " xarray.plot.plot(drill_cube.isel(x=[0],y=[0]).interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", " datasetID = datasetID + 1\n", " \n", "\n", " if geo_json['geometry']['type'] == 'Polygon':\n", "\n", " geom = transform_from_wgs_poly(geo_json['geometry'],int(EPSG))\n", "\n", " mask = rasterio.features.geometry_mask([geom for geoms in [geom]],\n", " out_shape=ds.geobox.shape,\n", " transform=ds.geobox.affine,\n", " all_touched=False,\n", " invert=True)\n", " \n", " dataMasked = target_dataset.where(mask)\n", " #data_masked_max = dataMasked.resample(time='1M').max()\n", " \n", " orig = dse.where(mask)\n", " #orig_max = orig.resample(time='1M').median()\n", " \n", " no_data_mask = orig.red.where(orig.red==-9999)\n", " no_data = no_data_mask.count(dim = ['x','y'])\n", " \n", " ## Prepare data for ploting extents\n", " cleaned_mask = no_data_mask.combine_first(dataMasked)\n", " \n", " nan_mask = mask.astype('float')\n", " nan_mask[~mask] = np.nan\n", " \n", " nans = np.isnan(nan_mask)\n", " nancols = np.all(nans, axis=0) \n", " nanrows = np.all(nans, axis=1) \n", " \n", " firstcol = nancols.argmin() \n", " firstrow = nanrows.argmin() \n", "\n", " lastcol = len(nancols) - nancols[::-1].argmin() # 8, last index where not NAN\n", " lastrow = len(nanrows) - nanrows[::-1].argmin() # 10\n", "\n", " # Area in Hectares\n", " no_data = (no_data * 625) / 10000\n", " total_area = (mask.sum() * 625) /10000\n", " \n", " data_masked = dataMasked.sum(dim = ['x','y'])\n", " \n", " other_masked = total_area - (data_masked + no_data)\n", "\n", " label = str(indices_buttons.value) + ' polygon ' + str(datasetID)\n", " labels.append(label)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", "\n", "\n", " geoJSONs.append(geo_json)\n", "\n", " fig = plt.figure(num=collect + 1, figsize=(10, 3.5))\n", " time_from = np.ravel(data_masked.time)\n", " \n", " plt.ylabel('Area (Ha)')\n", " color_scheme = '#1f77b4','#2ca02c', '#ff0000' ,'#d62728', \n", " plt.set_title = label\n", " #plt.set_yscale('log')\n", " plt.stackplot(time_from,\n", " data_masked.wofs.values, \n", " other_masked.wofs.values, \n", " no_data.values, \n", " colors = color_scheme,\n", " labels=['Water','Dry','No Data'])\n", " plt.legend(loc='upper left')\n", " fig.canvas.mpl_connect('button_press_event', onclick)\n", " collect = collect + 1 \n", " \n", " output_pandas = pd.concat([output_pandas,pd.DataFrame(data_masked.wofs.values,index = time_from, columns = [label])],axis = 1)\n", " #output_pandas.dropna(thresh=datasetID,inplace=True)\n", "\n", " out.clear_output()\n", " with out:\n", " ##Calc Stats.\n", " display(\"Mean Amount of Area (Ha) Classified as water over time period\")\n", " display(output_pandas.mean())\n", " display(\"Max Amount of Area (Ha) Classified as water over time period\")\n", " display(output_pandas.max())\n", " display(\"Min Amount of Area (Ha) Classified as water over time period\")\n", " display(output_pandas.min())\n", " \n", "\n", "draw_control3.on_draw(plot_nDVI)\n", "m3.add_control(draw_control3)\n", "\n", "label_html = widgets.HTML(\n", " value=\"Date selected from stack plot is \",\n", " layout=widgets.Layout(flex='1 1 auto', width='auto'))\n", "\n", "#box_layout = Layout(display='flex',\n", "# flex_flow='column',\n", "# align_items='stretch',\n", "# border='solid',\n", "# width='50%')\n", "\n", "button_plot = widgets.Button(description='Plot Spatial Extent', \n", " disabled=False,\n", " button_style='', # 'success', 'info', 'warning', 'danger' or ''\n", " tooltip='Click me', layout=widgets.Layout(flex='1 1 auto', width='auto'))\n", "\n", "button_plot.on_click(display_extent)\n", "\n", "\n", "plot_module = widgets.HBox([label_html, button_plot])\n", "\n", "pixel_drill = widgets.VBox([info,indices_buttons,m3,out,plot_module])\n", "\n", "pixel_drill\n", " " ] } ], "source": [ "%macro water_stackplot_aws_landsat8 8" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stored 'water_stackplot_aws_landsat8' (Macro)\n" ] } ], "source": [ "%store water_stackplot_aws_landsat8 " ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n", "No Case Study Selected\n", "No index selected\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "e103496332bc48db80a1ed84eb1ffa7c", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Output(layout=Layout(border='1px solid black')), ToggleButtons(description='Product:', options=…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "##Landsat 8 AWS WOFS Time series!!!\n", "%pylab notebook\n", "from importlib import reload # Python 3.4+ only.\n", "import sys\n", "import os\n", "import osr \n", "import ogr\n", "import ipywidgets as widgets\n", "import pandas as pd\n", "import rasterio.features\n", "import ipywidgets as widgets\n", "import xarray\n", "from utils import BandIndices\n", "\n", "#sys.path.append(os.path.expanduser('dea-notebooks/Scripts'))\n", "\n", "from utils import Water_Classifier\n", "Water_Classifier = reload(Water_Classifier)\n", "from utils.utils import (transform_from_wgs,\n", " transform_from_wgs_poly)\n", "import json\n", "from datacube import Datacube\n", "from datacube.storage import masking # Import masking capabilities\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "if 'case_study' not in locals():\n", " global case_study\n", " case_study = './configIndex.txt'\n", "\n", "\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "\n", "global output_pandas\n", "output_pandas= pd.DataFrame() \n", "\n", "global geoJSONs\n", "geoJSONs = []\n", "\n", "global datasetID\n", "datasetID = 0\n", "\n", "###Load EPSG\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", "\n", "##Define output widget to display table\n", "out = widgets.Output(layout={'border': '1px solid black'})\n", "info = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "def load_into_cube_aws_landsat8(measurements = None, cloud_mask = False):\n", "\n", " def lat_lon_to_epsg(lat_max, lon_min):\n", " return str(int(32700 - round((45 + lat_max) / 90, 0) * 100 + round((183 + lon_min) / 6, 0)))\n", "\n", " \n", " data = load_config(case_study)\n", "\n", " \n", " lats= data['lat']\n", " longs= data['lon']\n", " global EPSG\n", " EPSG = lat_lon_to_epsg(lats[1],longs[0])\n", " \n", " ds = dc.load(product='ls8_level1_usgs', \n", " measurements = ('red','green','blue','nir','swir1','swir2','quality'),\n", " group_by='solar_day',\n", " y=tuple(data['lat']),\n", " x=tuple(data['lon']),\n", " time=tuple(data['time']),\n", " output_crs='epsg:' + EPSG,\n", " resolution = (-25,25))\n", "\n", " #filter for clouds\n", " clean_pixel_mask = masking.make_mask(\n", " ds.quality,\n", " cloud=False,\n", " radiometric_saturation='none',\n", " terrain_occlusion = False)\n", " masked_cloud = ds.where(clean_pixel_mask)\n", "\n", "\n", " return masked_cloud\n", "\n", "#### Get geometry from config file draw polygon extent.\n", "config =load_config(case_study)\n", "\n", "geoJSON_Extent = config['geoJSON']\n", "#m3.add_layer(GeoJSON(data=draw_control2.last_draw))\n", "\n", "info_last_geom = ogr.CreateGeometryFromJson(str(geoJSON_Extent['geometry']))\n", "if info_last_geom:\n", " zoom = 13\n", " center = [info_last_geom.Centroid().GetY(),info_last_geom.Centroid().GetX()]\n", " m3 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "else:\n", " m3 = Map(center=center, zoom=zoom)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "\n", "dc = Datacube()\n", "\n", "indicesa = ['NDWI', 'WOFS']\n", "\n", "\n", "indices_buttons = widgets.ToggleButtons(\n", " options=indicesa,\n", " description='Product:',\n", " disabled=False,\n", " button_style='' \n", " #tooltips=indices_description,\n", " )\n", "\n", "\n", "if 'loading' in locals():\n", " ds = None\n", " print ('No Case Study Selected')\n", " target_dataset,desc = None, 'No index selected'\n", "else:\n", " geoJSON_Extent['properties']['style']['fillOpacity'] = 0\n", " geoJSON_Extent['properties']['style']['color'] = 'red'\n", " \n", " if 'loaded' in locals():\n", " if loaded != config['load_id']:\n", " print ('Retrieving New Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", " else:\n", " print ('Retrieving Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999) \n", " # Load data the first time.\n", "# indices_func(True)\n", " \n", "\n", " m3.add_layer(GeoJSON(data=geoJSON_Extent))\n", "\n", " \n", "def indices_func(change):\n", " global target_dataset\n", " if indices_buttons.value == 'NDWI':\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", "\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", "\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", "\n", " else:\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999)\n", " target_dataset,desc = Water_Classifier.water_classifier(dse)\n", "\n", " info.clear_output()\n", " with info:\n", " print(desc)\n", "\n", "\n", "with info:\n", " print(desc)\n", " \n", "indices_buttons.observe(indices_func, names = 'value')\n", "\n", "\n", "dataP = load_config(case_study)\n", "\n", "#EPSG = int(dataP['output_crs'].split(':')[1])\n", "#EPSG = 3577\n", "labels = []\n", "def plot_nDVI(self, action, geo_json):\n", " global output_pandas\n", " global datasetID\n", "\n", " #print (geo_json)\n", " if geo_json['geometry']['type'] == 'Point':\n", "\n", " pixel_drill_long,pixel_drill_lat = geo_json['geometry']['coordinates']\n", "\n", " pixel_drill_x, pixel_drill_y = transform_from_wgs(pixel_drill_long,pixel_drill_lat,int(EPSG))\n", " drill_cube = target_dataset.sel(x=[pixel_drill_x,pixel_drill_x + 25], y=[pixel_drill_y, pixel_drill_y + 25],\n", " method = 'nearest')\n", "\n", " time = np.ravel(drill_cube.isel(x=[0],y=[0]).time)\n", " values = np.ravel(drill_cube.isel(x=[0],y=[0]).values)\n", "\n", " #drill_cube.isel(x=[0],y=[0]).plot()\n", " label = str(indices_buttons.value) + ' point ' + str(datasetID)\n", " labels.append(label)\n", " \n", " xarray.plot.plot(drill_cube.isel(x=[0],y=[0]).interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", " datasetID = datasetID + 1\n", " \n", "\n", " if geo_json['geometry']['type'] == 'Polygon':\n", "\n", " geom = transform_from_wgs_poly(geo_json['geometry'],int(EPSG))\n", "\n", " mask = rasterio.features.geometry_mask([geom for geoms in [geom]],\n", " out_shape=ds.geobox.shape,\n", " transform=ds.geobox.affine,\n", " all_touched=False,\n", " invert=True)\n", "\n", " dataMasked = target_dataset.where(mask)\n", " \n", " #cleaned_dataset_month = dataMasked.resample(time='1M').median()\n", " orig = dse.where(mask)\n", " #no_data = orig.red.where(orig.red==-9999).count(dim = ['x','y'])\n", " \n", " cleaned_dataset_month = dataMasked.where(orig.red!=-9999)\n", " \n", " nan_mask = mask.astype('float')\n", " nan_mask[~mask] = np.nan\n", " \n", " nans = np.isnan(nan_mask)\n", " nancols = np.all(nans, axis=0) \n", " nanrows = np.all(nans, axis=1) \n", " \n", " firstcol = nancols.argmin() \n", " firstrow = nanrows.argmin() \n", "\n", " lastcol = len(nancols) - nancols[::-1].argmin() # 8, last index where not NAN\n", " lastrow = len(nanrows) - nanrows[::-1].argmin() # 10\n", " \n", " #no_data = no_data * 625\n", " total_area = mask.sum() * 625\n", " \n", "# data_masked = dataMasked.sum(dim = 'time')\n", " \n", " min_extent_time_ind = cleaned_dataset_month.wofs.where(cleaned_dataset_month.wofs == 0)\\\n", " .count(dim=['x','y']).argmax().values\n", " min_extent_acq = cleaned_dataset_month.wofs.isel(time=min_extent_time_ind)\n", " time_min = cleaned_dataset_month.time.isel(time=min_extent_time_ind)\n", " time_min_string = str(time_min.values)[0:10]\n", " min_extent_acq = min_extent_acq[firstrow:lastrow,firstcol:lastcol]\n", " \n", " # Find the acquisition with the most water.\n", " max_extent_time_ind = cleaned_dataset_month.wofs.where(cleaned_dataset_month.wofs > 0)\\\n", " .count(dim=['x','y']).argmax().values\n", " max_extent_acq = cleaned_dataset_month.wofs.isel(time=max_extent_time_ind)\n", " max_extent_acq = max_extent_acq[firstrow:lastrow,firstcol:lastcol]\n", " time_max = cleaned_dataset_month.time.isel(time=max_extent_time_ind)\n", " time_max_string = str(time_max.values)[0:10]\n", " \n", " #other_masked = total_area - (data_masked + no_data)\n", " \n", " #data_masked_mean.no_data = data_masked_mean - \n", " #data_masked_mean.plot()\n", "\n", " label = str(indices_buttons.value) + ' polygon ' + str(datasetID)\n", " labels.append(label)\n", "\n", " #xarray.plot.plot(data_masked_mean.interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " #plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", "\n", " datasetID = datasetID + 1\n", "\n", " geoJSONs.append(geo_json)\n", "\n", " \n", "# if first_plot == True:\n", " global f,axarr\n", "\n", " f, axarr = plt.subplots(2, figsize=(10, 10))\n", " #f.suptitle('Fractional Cover Landsat 8')\n", "\n", " axarr[0].set_title( 'Minimum Water Extent ' + time_min_string)\n", " axarr[1].set_title( 'Maximum Water Extent ' + time_max_string)\n", " #axarr[2].set_title( 'Mean Water Extent')\n", " #axarr[3].set_title( 'Class Extent Change')\n", " from matplotlib.colors import from_levels_and_colors\n", " cmap, norm = from_levels_and_colors([0,0.04,2],['red','blue'])\n", "\n", " axarr[0].imshow(min_extent_acq,cmap=cmap, norm=norm)\n", " axarr[1].imshow(max_extent_acq, cmap=cmap, norm=norm)\n", " \n", " #output_pandas = pd.concat([output_pandas,pd.DataFrame(data_masked.wofs.values,index = time_from, columns = [label])],axis = 1)\n", " #output_pandas.dropna(thresh=datasetID,inplace=True)\n", " #out.clear_output()\n", "# with out:\n", "# display(output_pandas)\n", "\n", "draw_control3.on_draw(plot_nDVI)\n", "m3.add_control(draw_control3)\n", "\n", "pixel_drill = widgets.VBox([info,indices_buttons,m3,out])\n", "\n", "pixel_drill" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Macro `water_max_min_aws_landsat8` created. To execute, type its name (without quotes).\n", "=== Macro contents: ===\n", "##Landsat 8 AWS WOFS Time series!!!\n", "get_ipython().run_line_magic('pylab', 'notebook')\n", "from importlib import reload # Python 3.4+ only.\n", "import sys\n", "import os\n", "import osr \n", "import ogr\n", "import ipywidgets as widgets\n", "import pandas as pd\n", "import rasterio.features\n", "import ipywidgets as widgets\n", "import xarray\n", "from utils import BandIndices\n", "\n", "#sys.path.append(os.path.expanduser('dea-notebooks/Scripts'))\n", "\n", "from utils import Water_Classifier\n", "Water_Classifier = reload(Water_Classifier)\n", "from utils.utils import (transform_from_wgs,\n", " transform_from_wgs_poly)\n", "import json\n", "from datacube import Datacube\n", "from datacube.storage import masking # Import masking capabilities\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "if 'case_study' not in locals():\n", " global case_study\n", " case_study = './configIndex.txt'\n", "\n", "\n", "from ipyleaflet import (\n", " Map,\n", " Marker,\n", " TileLayer, ImageOverlay,\n", " Polyline, Polygon, Rectangle, Circle, CircleMarker,\n", " GeoJSON,\n", " DrawControl,\n", " basemaps\n", ")\n", "\n", "global output_pandas\n", "output_pandas= pd.DataFrame() \n", "\n", "global geoJSONs\n", "geoJSONs = []\n", "\n", "global datasetID\n", "datasetID = 0\n", "\n", "###Load EPSG\n", "def load_config(filename):\n", " with open(filename, 'r') as f:\n", " data = json.load(f)\n", " return data\n", "\n", "##Define output widget to display table\n", "out = widgets.Output(layout={'border': '1px solid black'})\n", "info = widgets.Output(layout={'border': '1px solid black'})\n", "\n", "\n", "def load_into_cube_aws_landsat8(measurements = None, cloud_mask = False):\n", "\n", " def lat_lon_to_epsg(lat_max, lon_min):\n", " return str(int(32700 - round((45 + lat_max) / 90, 0) * 100 + round((183 + lon_min) / 6, 0)))\n", "\n", " \n", " data = load_config(case_study)\n", "\n", " \n", " lats= data['lat']\n", " longs= data['lon']\n", " global EPSG\n", " EPSG = lat_lon_to_epsg(lats[1],longs[0])\n", " \n", " ds = dc.load(product='ls8_level1_usgs', \n", " measurements = ('red','green','blue','nir','swir1','swir2','quality'),\n", " group_by='solar_day',\n", " y=tuple(data['lat']),\n", " x=tuple(data['lon']),\n", " time=tuple(data['time']),\n", " output_crs='epsg:' + EPSG,\n", " resolution = (-25,25))\n", "\n", " #filter for clouds\n", " clean_pixel_mask = masking.make_mask(\n", " ds.quality,\n", " cloud=False,\n", " radiometric_saturation='none',\n", " terrain_occlusion = False)\n", " masked_cloud = ds.where(clean_pixel_mask)\n", "\n", "\n", " return masked_cloud\n", "\n", "#### Get geometry from config file draw polygon extent.\n", "config =load_config(case_study)\n", "\n", "geoJSON_Extent = config['geoJSON']\n", "#m3.add_layer(GeoJSON(data=draw_control2.last_draw))\n", "\n", "info_last_geom = ogr.CreateGeometryFromJson(str(geoJSON_Extent['geometry']))\n", "if info_last_geom:\n", " zoom = 13\n", " center = [info_last_geom.Centroid().GetY(),info_last_geom.Centroid().GetX()]\n", " m3 = Map(center=center, zoom=zoom, basemap=basemaps.Esri.WorldImagery)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "else:\n", " m3 = Map(center=center, zoom=zoom)\n", " draw_control3 = DrawControl(polygon = {\"shapeOptions\": {\"fillOpacity\": 0}})\n", "\n", "dc = Datacube()\n", "\n", "indicesa = ['NDWI', 'WOFS']\n", "\n", "\n", "indices_buttons = widgets.ToggleButtons(\n", " options=indicesa,\n", " description='Product:',\n", " disabled=False,\n", " button_style='' \n", " #tooltips=indices_description,\n", " )\n", "\n", "\n", "if 'loading' in locals():\n", " ds = None\n", " print ('No Case Study Selected')\n", " target_dataset,desc = None, 'No index selected'\n", "else:\n", " geoJSON_Extent['properties']['style']['fillOpacity'] = 0\n", " geoJSON_Extent['properties']['style']['color'] = 'red'\n", " \n", " if 'loaded' in locals():\n", " if loaded != config['load_id']:\n", " print ('Retrieving New Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", " else:\n", " print ('Retrieving Cube')\n", " ds = load_into_cube_aws_landsat8()\n", " loaded = config['load_id']\n", "\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999) \n", " # Load data the first time.\n", "# indices_func(True)\n", " \n", "\n", " m3.add_layer(GeoJSON(data=geoJSON_Extent))\n", "\n", " \n", "def indices_func(change):\n", " global target_dataset\n", " if indices_buttons.value == 'NDWI':\n", " water_dataset, desc = BandIndices.calculate_indices(ds,indices_buttons.value)\n", "\n", " target_dataset_1 = (water_dataset > 0.05).astype(np.float32)\n", "\n", " target_dataset = target_dataset_1.to_dataset(name='wofs')\n", "\n", " else:\n", " #dse = ds.resample(time='1M').median()\n", " dse = ds.fillna(-9999)\n", " target_dataset,desc = Water_Classifier.water_classifier(dse)\n", "\n", " info.clear_output()\n", " with info:\n", " print(desc)\n", "\n", "\n", "with info:\n", " print(desc)\n", " \n", "indices_buttons.observe(indices_func, names = 'value')\n", "\n", "\n", "dataP = load_config(case_study)\n", "\n", "#EPSG = int(dataP['output_crs'].split(':')[1])\n", "#EPSG = 3577\n", "labels = []\n", "def plot_nDVI(self, action, geo_json):\n", " global output_pandas\n", " global datasetID\n", "\n", " #print (geo_json)\n", " if geo_json['geometry']['type'] == 'Point':\n", "\n", " pixel_drill_long,pixel_drill_lat = geo_json['geometry']['coordinates']\n", "\n", " pixel_drill_x, pixel_drill_y = transform_from_wgs(pixel_drill_long,pixel_drill_lat,int(EPSG))\n", " drill_cube = target_dataset.sel(x=[pixel_drill_x,pixel_drill_x + 25], y=[pixel_drill_y, pixel_drill_y + 25],\n", " method = 'nearest')\n", "\n", " time = np.ravel(drill_cube.isel(x=[0],y=[0]).time)\n", " values = np.ravel(drill_cube.isel(x=[0],y=[0]).values)\n", "\n", " #drill_cube.isel(x=[0],y=[0]).plot()\n", " label = str(indices_buttons.value) + ' point ' + str(datasetID)\n", " labels.append(label)\n", " \n", " xarray.plot.plot(drill_cube.isel(x=[0],y=[0]).interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", " datasetID = datasetID + 1\n", " \n", "\n", " if geo_json['geometry']['type'] == 'Polygon':\n", "\n", " geom = transform_from_wgs_poly(geo_json['geometry'],int(EPSG))\n", "\n", " mask = rasterio.features.geometry_mask([geom for geoms in [geom]],\n", " out_shape=ds.geobox.shape,\n", " transform=ds.geobox.affine,\n", " all_touched=False,\n", " invert=True)\n", "\n", " dataMasked = target_dataset.where(mask)\n", " \n", " #cleaned_dataset_month = dataMasked.resample(time='1M').median()\n", " orig = dse.where(mask)\n", " #no_data = orig.red.where(orig.red==-9999).count(dim = ['x','y'])\n", " \n", " cleaned_dataset_month = dataMasked.where(orig.red!=-9999)\n", " \n", " nan_mask = mask.astype('float')\n", " nan_mask[~mask] = np.nan\n", " \n", " nans = np.isnan(nan_mask)\n", " nancols = np.all(nans, axis=0) \n", " nanrows = np.all(nans, axis=1) \n", " \n", " firstcol = nancols.argmin() \n", " firstrow = nanrows.argmin() \n", "\n", " lastcol = len(nancols) - nancols[::-1].argmin() # 8, last index where not NAN\n", " lastrow = len(nanrows) - nanrows[::-1].argmin() # 10\n", " \n", " #no_data = no_data * 625\n", " total_area = mask.sum() * 625\n", " \n", "# data_masked = dataMasked.sum(dim = 'time')\n", " \n", " min_extent_time_ind = cleaned_dataset_month.wofs.where(cleaned_dataset_month.wofs == 0) .count(dim=['x','y']).argmax().values\n", " min_extent_acq = cleaned_dataset_month.wofs.isel(time=min_extent_time_ind)\n", " time_min = cleaned_dataset_month.time.isel(time=min_extent_time_ind)\n", " time_min_string = str(time_min.values)[0:10]\n", " min_extent_acq = min_extent_acq[firstrow:lastrow,firstcol:lastcol]\n", " \n", " # Find the acquisition with the most water.\n", " max_extent_time_ind = cleaned_dataset_month.wofs.where(cleaned_dataset_month.wofs > 0) .count(dim=['x','y']).argmax().values\n", " max_extent_acq = cleaned_dataset_month.wofs.isel(time=max_extent_time_ind)\n", " max_extent_acq = max_extent_acq[firstrow:lastrow,firstcol:lastcol]\n", " time_max = cleaned_dataset_month.time.isel(time=max_extent_time_ind)\n", " time_max_string = str(time_max.values)[0:10]\n", " \n", " #other_masked = total_area - (data_masked + no_data)\n", " \n", " #data_masked_mean.no_data = data_masked_mean - \n", " #data_masked_mean.plot()\n", "\n", " label = str(indices_buttons.value) + ' polygon ' + str(datasetID)\n", " labels.append(label)\n", "\n", " #xarray.plot.plot(data_masked_mean.interpolate_na(dim = 'time', method = 'nearest'), marker='*')\n", " #plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=4, mode=\"expand\", labels = labels)\n", "\n", " geo_json['properties']['datasetID'] = datasetID\n", "\n", " datasetID = datasetID + 1\n", "\n", " geoJSONs.append(geo_json)\n", "\n", " \n", "# if first_plot == True:\n", " global f,axarr\n", "\n", " f, axarr = plt.subplots(2, figsize=(10, 10))\n", " #f.suptitle('Fractional Cover Landsat 8')\n", "\n", " axarr[0].set_title( 'Minimum Water Extent ' + time_min_string)\n", " axarr[1].set_title( 'Maximum Water Extent ' + time_max_string)\n", " #axarr[2].set_title( 'Mean Water Extent')\n", " #axarr[3].set_title( 'Class Extent Change')\n", " from matplotlib.colors import from_levels_and_colors\n", " cmap, norm = from_levels_and_colors([0,0.04,2],['red','blue'])\n", "\n", " axarr[0].imshow(min_extent_acq,cmap=cmap, norm=norm)\n", " axarr[1].imshow(max_extent_acq, cmap=cmap, norm=norm)\n", " \n", " #output_pandas = pd.concat([output_pandas,pd.DataFrame(data_masked.wofs.values,index = time_from, columns = [label])],axis = 1)\n", " #output_pandas.dropna(thresh=datasetID,inplace=True)\n", " #out.clear_output()\n", "# with out:\n", "# display(output_pandas)\n", "\n", "draw_control3.on_draw(plot_nDVI)\n", "m3.add_control(draw_control3)\n", "\n", "pixel_drill = widgets.VBox([info,indices_buttons,m3,out])\n", "\n", "pixel_drill\n", " " ] } ], "source": [ "%macro water_max_min_aws_landsat8 11" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stored 'water_max_min_aws_landsat8' (Macro)\n" ] } ], "source": [ "%store water_max_min_aws_landsat8 " ] } ], "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.6.7" } }, "nbformat": 4, "nbformat_minor": 2 }