<jupyter_start><jupyter_text>A quick and dirty LandSAT validation testIn this notebook we extract a subset of Australian tree cover from LandSAT data, and do a quick and dirty validation exercise with coincident LiDAR.The point of this exercise is to test and demonstrate data service models, *not* to perform an accurate validation.<jupyter_code># import python things - we will need:

#OS interaction (maybe)
import os

#postgres 
import psycopg2 as ppg

#numpy
import numpy as np

#osgeo libraries
from osgeo import ogr

#possibly - shapely for tinkering with geometries
import shapely as sp
from shapely.geometry import Point,Polygon,asShape
from shapely.wkt import loads as wkt_loads
from shapely import speedups

#maybe cartopy, but I'm not sure
#import cartopy as cp
#import cartopy.crs as ccrs

#geopandas - very neat for pulling up postGIS geometries
import geopandas as gp

#matplotlib of course!
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib.collections import PatchCollection
from mpl_toolkits.basemap import Basemap

#set up for 3d plots
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pylab as pylab

#fiona and descartes
import fiona
from descartes import PolygonPatch

#OWSlib, so we can get web coverages etc.
from owslib.wms import WebMapService
from owslib.wcs import WebCoverageService


import numpy as np
from netCDF4 import Dataset
import pandas as pd

z1 = np.greater(2,1)
z1<jupyter_output><empty_output><jupyter_text>1. Define a region and select some LandSAT data. We know our LIDAR data only cover a small region:...and were collected in 2013 - so we'll grab coincident LandSAT data from the same year and the same region using a web service request<jupyter_code>f = Dataset('http://dapds00.nci.org.au/thredds/dodsC/ub8/au/treecover/ANUWALD.TreeCover.25m.2015.nc')
#don't need to pull any actual data down yet - lets see what is in this file.
for key in f.variables.keys():
    print(f.variables[key])
lon = f.variables['longitude'][:]
lat = f.variables['latitude'][:]
##so lets find a Canberra-ish region
region = [134,136, -38, -36]
lon_si = np.squeeze(np.where((lon >= region[0]) & (lon <= region[1])))
lat_si = np.squeeze(np.where((lat >= region[2]) & (lat <= region[3])))
lat[lat_si[0]]
(min(lon_si),max(lon_si))
lon_si.shape
##and harvest the region from the netCDF file
# can we directly use lon and lat on the GVC array? hmm
# No, but we don't have to get the array before indexing into it.. we just get the subset!
#treecover = f.variables['TreeCover'][min(lon_si[0]): max(lon_si[0]), min(lat_si[0]): max(lat_si[0])]

treecover = f.variables['TreeCover'][0,110000:129999,92000:111999]
tree_cover = treecover[0]
tree_cover = np.ma.filled(tree_cover, -999)
tree_cover.shape
tree_cover = np.ma.filled(treecover)
tree_cover
from cartopy import config
import cartopy.crs as ccrs
# Plot. 
fig = plt.figure(figsize=(10, 10))

plt.imshow(tree_cover)
this = np.where((tree_cover > 0) & (tree_cover <= 1))
np.shape(tree_cover)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09381107491856677, 'fraction_numerical': 0.025732899022801303, 'mean_word_length': 5.305954825462012, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '266632', 'license': 'Apache-2.0', 'n_tokens_mistral': 1114, 'n_tokens_neox': 1048, 'n_words': 363, 'path': 'pgpointcloud/lidar_landsat.ipynb', 'repo_name': 'adamsteer/nci-notebooks', 'size': '3070'}