<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 How to make a 3D web visualisation without a single line of code
In this notebook we use QGIS to create a shareable terrain model with a data overlay, which can be shared on a web server, without typing a single line of code. Let's go!
To inetract with the maps we make below, check this notebook out in the ipython notebook viewer
Step1:  and there - your first 3D interactive map, made with no coding and using web data services!
...but it isn't very useful or informative. How can we fix that?
5. Add some more data
We might be interested in land cover attributes for our region, so let's get them! How about photosynthetic vegetation for 2015?
http
Step2:  Now we have an elevation map coloured by green vegetation! But that's only a pretty picture.
6. Now lets do some analysis!
We will try to show
Step3:  Now, let's add some complexity
We have some elevation data, we have some cadastral data, we have some data about photosynthetic vegetation cover. We can do our quick visualisation of whether block hilliness (as determined by SRTM height standard deviation for each block) is related to photosynthetic vegetation cover (as determined by the median of vegetation cover inside each block). We can set this up as follows
<ASSISTANT_TASK:>
Python Code:
###ignore this block of code - it is required only to show the map in iPython - you won't need it!
from IPython.core.display import display, HTML
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/ACT_elevs_test_1.html"></iframe>'))
Explanation: How to make a 3D web visualisation without a single line of code
In this notebook we use QGIS to create a shareable terrain model with a data overlay, which can be shared on a web server, without typing a single line of code. Let's go!
To inetract with the maps we make below, check this notebook out in the ipython notebook viewer:
http://nbviewer.jupyter.org/github/adamsteer/nci-notebooks/blob/master/How%20to%20make%203D%20visualisations%20from%20NCI%20data%20without%20any%20coding%20%28nearly%29.ipynb
1. Decide on some datasets
Let's use Canberra. We need:
* a terrain model
* some data - what about... false colour LandSAT imagery? or NVDI? or how about vegetation types?
* some more data - how about OpenStreetMap buildings?
2. Get the topography data and import into QGIS
In this approach we download all our data from NCI to our machine (or could use it in the VDI), because we need to modify bands - which QGIS is not happy to do for web services (WMS or WCS).
For topography we need a terrain model in a raster format, e.g. GeoTIFF, which covers Canberra. We can head to the NCI Elevation collection here:
http://dapds00.nci.org.au/thredds/catalog/rr1/Elevation/catalog.html
...and look for SRTM 1 second elevation - good enough for this job. If you are happy with ESRI grids, navigate to the tile collection here:
http://dapds00.nci.org.au/thredds/catalog/rr1/Elevation/1secSRTM_DEMs_v1.0/DEM/Tiles/catalog.html
Data are organised in folders named by longitide and latitude of the south west (bottom left) corner. For Canberra we need 149, -36.
...and for now, this is probably the best method. We could make an OpenDAP or WCS request for a subset, but that would be coding! The pull of the dark side is strong - so here is a link that gets a GeoTIFF from the SRTM tile, using WCS:
http://dapds00.nci.org.au/thredds/wcs/rr1/Elevation/NetCDF/1secSRTM_DEMs_v1.0/DEM/Elevation_1secSRTM_DEMs_v1.0_DEM_Tiles_e149s36dem1_0.nc?service=WCS&version=1.0.0&request=GetCoverage&Coverage=elevation&bbox=149.0,-36,149.9,-35&format=GeoTIFF
Now, import the resulting GeoTIFF into QGIS as a raster data source:
3. Install a QGIS plugin and use it.
This is seriously not coding - just head to the "plugins" menu, click "manage and install plugins", and find the Qgis2threejs plugin. Install it:
4. Set up your first interactive map
Now zoom into the DEM so that it covers the entire display window, then head to the 'web' menu. Choose 'qgis2threejs'. In the resulting dialog box, click 'world' in the left pane, and find 'vertical exaggeration'. Set it to 10.
In 'DEM', the one active layer (your GeoTIFF) should be preselected. Click 'run', and you should get a web browser opening with a 3D model inside!
The expected output is shown below - try viewing the notebook here:
http://nbviewer.jupyter.org/github/adamsteer/nci-notebooks/blob/master/How%20to%20make%203D%20visualisations%20from%20NCI%20data%20without%20any%20coding%20%28nearly%29.ipynb
...you should be able to move the map, zoom in and out, and generally inspect a terrain model.
End of explanation
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/act_elevs_plus_greenveg.html"></iframe>'))
Explanation: and there - your first 3D interactive map, made with no coding and using web data services!
...but it isn't very useful or informative. How can we fix that?
5. Add some more data
We might be interested in land cover attributes for our region, so let's get them! How about photosynthetic vegetation for 2015?
http://dapds00.nci.org.au/thredds/catalog/ub8/au/FractCov/PV/catalog.html?dataset=ub8-au/FractCov/PV/FractCover.V3_0_1.2015.aust.005.PV.nc
Using WCS again - click on the WCS link, and look for a <name> tag - it says 'PV'. This is the coverage we need to get. So we form a WCS request like this:
the dataset: http://dapds00.nci.org.au/thredds/wcs/ub8/au/FractCov/PV/FractCover.V3_0_1.2015.aust.005.PV.nc
the service: service=WCS
the service version: version=1.0.0
the thing we want to do (get a coverage): request=GetCoverage
the coverage (or layer) we want to get: Coverage=PV
the boundary of the layer we want: bbox=149.0,-36,149.9,-35
the format we want to get our coverage as: format=GeoTIFF
...so we put a question mark after the dataset name, then add the rest of the labels describing the thing we want afterward, in any order, separated by ampersands:
http://dapds00.nci.org.au/thredds/wcs/ub8/au/FractCov/PV/FractCover.V3_0_1.2015.aust.005.PV.nc?service=WCS&version=1.0.0&request=GetCoverage&Coverage=PV&bbox=149.0,-36,149.9,-35&format=GeoTIFF
(woah. There was a glitch in the matrix - if we didn't write out that URL, you would have had to code just now).
Add the resulting GeoTIFF to QGIS as a raster data source, just like the DEM! Once you have a style you're happy with, soom to the desired extent and use the qgis2threejs plugin to make a new map.
ACT DEM coloured by green vegetation content
Click on the map to inspect features, click and drag to move, scroll to zoom.
End of explanation
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/act_block_hilliness_proxy.html"></iframe>'))
Explanation: Now we have an elevation map coloured by green vegetation! But that's only a pretty picture.
6. Now lets do some analysis!
We will try to show:
Standard deviaton of elevation of blocks as a proxy for hilliness, plotted as a volume on the elevation map
Sum of fractional PV cover (maybe bare ground, or tree cover?) for each block, also as a volume on the map
Interactive layer selection using three.js maps
...and no code, Only button clicking.
Back to QGIS. We need another plugin - zonal statistics! We know how to install plugins, so get to it. When ready, we'll make a hilliness proxy first.
Next we will make a vegetation cover proxy
..and finally, can we come up with a metric of vegetation cover as a function of block hilliness? A question here might be 'do land owners tend to clear flatter blocks more than hilly blocks?'. Can we answer it using web services data, QGIS, some clicks - and then make a very pretty, and interactive map?
On to the first question - how do we get ACT block data? Head to ACTMAPi and find ACT blocks. Here's a shortcut:
http://actmapi.actgov.opendata.arcgis.com/datasets/afa1d909a0ae427cb9c1963e0d2e80ca_4
Find the 'API' menu box, click the down arrow and then copy the URL inside the JSON option:
http://actmapi.actgov.opendata.arcgis.com/datasets/afa1d909a0ae427cb9c1963e0d2e80ca_4.geojson
In QGIS, add a new vector layer. Choose 'service' from the options, and paste the URL into the appropriate box:
Wait a while! Now because QGIS you'll need to save the result as a shapefile before we can do anything with it. Save the layer, close the GeoJSON layer and reopen your new ACT blocks shapefile.
...it's probably simpler to just download and open the shapefile - but now we know how to add vector layers from a web service> either way, the result is something like this:
We don't need to make the blocks pretty yet, We'll do something!
Block hilliness
Using the SRTM DEM, let's make a proxy of block hilliness. Check that you have the QGIS zonal statistics plugin (Raster -> zonal statistics). If not, install it the same way you installed qgis2three.js. Once you have it, open the plugin. Choose the DEM as the raster layer, and use band 1. Then choose your ACT blocks vector layer. In the statistics to calculate, pick an appropriate set - but include standard deviation - this is our roughness proxy.
Run the plugin, then open the properties box of the ACT blocks layer and colour your blocks by standard deviation. The ZOnal statistics plugin has looped through all the polygons in the blocks layer and computed descriptive statistics of the underlying DEM. And there you have it - ACT blocks coloured by the standard deviation of the elevation they contain, as a proxy for hilliness. 
Now lets make a cool map! zoom in so that you have a region you're happy with occupying your map view, and open the qgis2threejs plugin.
In the left pane, under polygon option choose your ACT blocks layer. Options for styling it appear on the right - which I'll recreate again because I crashed QGIS again at a bad time (Save, fella, save!)
Result - 3D block hilliness visualised as height of a block polygon
What we see here are cadastral blocks, with our 'hilliness proxy' displayed by colour and extruded column height. Darker blue, and taller columns are hillier! Click on the map to inspect features, click and drag to move, scroll to zoom.
End of explanation
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling ="no" src="./qgis2threejs/veg_mean_colours.html"></iframe>'))
Explanation: Now, let's add some complexity
We have some elevation data, we have some cadastral data, we have some data about photosynthetic vegetation cover. We can do our quick visualisation of whether block hilliness (as determined by SRTM height standard deviation for each block) is related to photosynthetic vegetation cover (as determined by the median of vegetation cover inside each block). We can set this up as follows:
classify and colour blocks by vegetation cover
visualise block hilliness as the height of an extruded column
In this scheme, if our hypothesis is that hillier blocks are less cleared, dark green blocks will visualise as taller columns. Lets test it out!
End of explanation
<END_TASK>