#!/usr/bin/env python
# coding: utf-8

# [En français](../tutorial_gdal_fr/)  
# 
# ![ECCC logo](https://eccc-msc.github.io/open-data/img_eccc-logo.png)  
# 
# [TOC](https://eccc-msc.github.io/open-data/readme_en/) > [Usage overview](https://eccc-msc.github.io/open-data/usage/readme_en/) > GDAL command line tutorial
# 
# # GDAL command line tutorial with weather data
# 
# ## Introduction
# 
# [MSC GeoMet](https://eccc-msc.github.io/open-data/msc-geomet/readme_en/) and [MSC Datamart](https://eccc-msc.github.io/open-data/msc-datamart/readme_en/) data can be manipulated via the command line using [GDAL](https://gdal.org/), a widely-known software library used to read and write raster and vector geospatial data. In the following examples, you'll use a GeoTIFF file retrieved using via a Web Coverage Service (WCS) request to MSC GeoMet. This tutorial will show you how to:
# * Display the GDAL version installed on your system
# * Save a WCS request output to your computer
# * List information/metadata related to the raster file
# * Reproject a raster file
# * Convert a GeoTIFF file to the NetCDF file format
# * Get the value for a specific point based on a location in longitude/latitude
# 
# There are various ways to install GDAL, please refer to https://gdal.org/ for more information.
# 
# To run the following GDAL command line examples, you need to have a basic knowledge of using the terminal command line. These examples work within a bash terminal. 
# 
# The [interactive version of this Jupyter Notebook is available](https://mybinder.org/v2/gh/ECCC-MSC/open-data/master?filepath=docs%2Fusage%2Ftutorial_gdal%2ftutorial_gdal_en.ipynb).
# 
# [![badge](https://img.shields.io/badge/Interactive%20version-binder-F5A252.svg?logo=data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFkAAABZCAMAAABi1XidAAAB8lBMVEX///9XmsrmZYH1olJXmsr1olJXmsrmZYH1olJXmsr1olJXmsrmZYH1olL1olJXmsr1olJXmsrmZYH1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olJXmsrmZYH1olL1olL0nFf1olJXmsrmZYH1olJXmsq8dZb1olJXmsrmZYH1olJXmspXmspXmsr1olL1olJXmsrmZYH1olJXmsr1olL1olJXmsrmZYH1olL1olLeaIVXmsrmZYH1olL1olL1olJXmsrmZYH1olLna31Xmsr1olJXmsr1olJXmsrmZYH1olLqoVr1olJXmsr1olJXmsrmZYH1olL1olKkfaPobXvviGabgadXmsqThKuofKHmZ4Dobnr1olJXmsr1olJXmspXmsr1olJXmsrfZ4TuhWn1olL1olJXmsqBi7X1olJXmspZmslbmMhbmsdemsVfl8ZgmsNim8Jpk8F0m7R4m7F5nLB6jbh7jbiDirOEibOGnKaMhq+PnaCVg6qWg6qegKaff6WhnpKofKGtnomxeZy3noG6dZi+n3vCcpPDcpPGn3bLb4/Mb47UbIrVa4rYoGjdaIbeaIXhoWHmZYHobXvpcHjqdHXreHLroVrsfG/uhGnuh2bwj2Hxk17yl1vzmljzm1j0nlX1olL3AJXWAAAAbXRSTlMAEBAQHx8gICAuLjAwMDw9PUBAQEpQUFBXV1hgYGBkcHBwcXl8gICAgoiIkJCQlJicnJ2goKCmqK+wsLC4usDAwMjP0NDQ1NbW3Nzg4ODi5+3v8PDw8/T09PX29vb39/f5+fr7+/z8/Pz9/v7+zczCxgAABC5JREFUeAHN1ul3k0UUBvCb1CTVpmpaitAGSLSpSuKCLWpbTKNJFGlcSMAFF63iUmRccNG6gLbuxkXU66JAUef/9LSpmXnyLr3T5AO/rzl5zj137p136BISy44fKJXuGN/d19PUfYeO67Znqtf2KH33Id1psXoFdW30sPZ1sMvs2D060AHqws4FHeJojLZqnw53cmfvg+XR8mC0OEjuxrXEkX5ydeVJLVIlV0e10PXk5k7dYeHu7Cj1j+49uKg7uLU61tGLw1lq27ugQYlclHC4bgv7VQ+TAyj5Zc/UjsPvs1sd5cWryWObtvWT2EPa4rtnWW3JkpjggEpbOsPr7F7EyNewtpBIslA7p43HCsnwooXTEc3UmPmCNn5lrqTJxy6nRmcavGZVt/3Da2pD5NHvsOHJCrdc1G2r3DITpU7yic7w/7Rxnjc0kt5GC4djiv2Sz3Fb2iEZg41/ddsFDoyuYrIkmFehz0HR2thPgQqMyQYb2OtB0WxsZ3BeG3+wpRb1vzl2UYBog8FfGhttFKjtAclnZYrRo9ryG9uG/FZQU4AEg8ZE9LjGMzTmqKXPLnlWVnIlQQTvxJf8ip7VgjZjyVPrjw1te5otM7RmP7xm+sK2Gv9I8Gi++BRbEkR9EBw8zRUcKxwp73xkaLiqQb+kGduJTNHG72zcW9LoJgqQxpP3/Tj//c3yB0tqzaml05/+orHLksVO+95kX7/7qgJvnjlrfr2Ggsyx0eoy9uPzN5SPd86aXggOsEKW2Prz7du3VID3/tzs/sSRs2w7ovVHKtjrX2pd7ZMlTxAYfBAL9jiDwfLkq55Tm7ifhMlTGPyCAs7RFRhn47JnlcB9RM5T97ASuZXIcVNuUDIndpDbdsfrqsOppeXl5Y+XVKdjFCTh+zGaVuj0d9zy05PPK3QzBamxdwtTCrzyg/2Rvf2EstUjordGwa/kx9mSJLr8mLLtCW8HHGJc2R5hS219IiF6PnTusOqcMl57gm0Z8kanKMAQg0qSyuZfn7zItsbGyO9QlnxY0eCuD1XL2ys/MsrQhltE7Ug0uFOzufJFE2PxBo/YAx8XPPdDwWN0MrDRYIZF0mSMKCNHgaIVFoBbNoLJ7tEQDKxGF0kcLQimojCZopv0OkNOyWCCg9XMVAi7ARJzQdM2QUh0gmBozjc3Skg6dSBRqDGYSUOu66Zg+I2fNZs/M3/f/Grl/XnyF1Gw3VKCez0PN5IUfFLqvgUN4C0qNqYs5YhPL+aVZYDE4IpUk57oSFnJm4FyCqqOE0jhY2SMyLFoo56zyo6becOS5UVDdj7Vih0zp+tcMhwRpBeLyqtIjlJKAIZSbI8SGSF3k0pA3mR5tHuwPFoa7N7reoq2bqCsAk1HqCu5uvI1n6JuRXI+S1Mco54YmYTwcn6Aeic+kssXi8XpXC4V3t7/ADuTNKaQJdScAAAAAElFTkSuQmCC)](https://mybinder.org/v2/gh/ECCC-MSC/open-data/master?filepath=docs%2Fusage%2Ftutorial_gdal%2ftutorial_gdal_en.ipynb)
# 
# ## Show GDAL version
# 
# GDAL is a suite of several command line tools. When you install GDAL you get all the different command line tools. The most basic tool is `gdalinfo` which can be used to retrieve information pertaining to your GDAL installation and display information about a raster file.

# In[1]:


get_ipython().run_cell_magic('bash', '', '\ngdalinfo --version\n')


# ## Save a WCS request output to disk
# 
# The OGC Web Coverage Service requests enable a client to retrieve coverage information from a raster file for a given area of interest. WCS requests are made over the internet (HTTPS) and give the user more flexibility when requesting information about the coverage of a layer compared with the more traditional way of downloading flat files. The Web Coverage Service allows for several different types of requests, each of which are described in further detail below. For more information on the WCS request parameters, please refer to the [MSC GeoMet WCS GetCoverage page](https://eccc-msc.github.io/open-data/msc-geomet/web-services_en/#wcs-getcoverage).
# 
# We are going to use a `curl` command to save the WCS request result on disk, the file will be named `CMC_glb_TMP.tif`. The result is a GeoTIFF file, showing temperature (°C) from a subset of MSC's Global Deterministic Prediction System (GDPS).

# In[4]:


get_ipython().run_cell_magic('bash', '', '\ncurl "https://geo.weather.gc.ca/geomet?SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=GDPS.ETA_TT&SUBSETTINGCRS=EPSG:4326&SUBSET=x(-120,-85)&SUBSET=y(48,66)&RESOLUTION=x(0.24)&RESOLUTION=y(0.24)&FORMAT=image/tiff" > CMC_glb_TMP.tif \n')


# ## Lists information about a raster file
# 
# The `gdalinfo` tool can be used to retrieve the downloaded raster file's metadata. The command's output will list some information on the file, such as:
# * File format
# * File size
# * Coordinate system
# * Metadata
# * Band information

# In[11]:


get_ipython().run_cell_magic('bash', '', '\ngdalinfo CMC_glb_TMP.tif\n')


# It is also possible to use `gdalinfo` to retrieve some basic statistics on the raster file, such as minimum and maximum value by adding the `-mm` option. Note that the resulting values are in °C.

# In[21]:


get_ipython().run_cell_magic('bash', '', '\ngdalinfo -mm CMC_glb_TMP.tif | grep Min/Max\n')


# Adding the `-proj4` option to `gdalinfo` will output the projection definition as a proj4 string:

# In[24]:


get_ipython().run_cell_magic('bash', '', '\ngdalinfo -proj4 CMC_glb_TMP.tif | grep PROJ.4 -A 1\n')


# ## Reproject a raster file
# 
# Using the `gdalwarp` command, we can reproject a raster file. Using MSC GeoMet and MSC Datamart data, you only need to provide an output projection definition corresponding to a EPSG code, or you can use a proj4 string.
# 
# The following example reprojects the GeoTIFF file into an EPSG:3857 projection. The output file is named `CMC_glb_TMP_epsg3857.tif`

# In[25]:


get_ipython().run_cell_magic('bash', '', '\ngdalwarp -t_srs EPSG:3857 CMC_glb_TMP.tif CMC_glb_TMP_epsg3857.tif\n')


# Then we can use `gdalinfo` to look at the coordinates and the proj4 string to ensure that the projection of `CMC_glb_TMP_epsg3857.tif` really is different from the original file.

# In[26]:


get_ipython().run_cell_magic('bash', '', "\ngdalinfo -proj4 epsg3857.tif | grep -E '(PROJ.4|Corner Coordinates:)' -A 5\n")


# ## Convert a GeoTIFF file to the NetCDF file format
# 
# Using `gdal_translate` command, we can convert a raster file from any supported format (`gdalinfo --formats`) into another raster file format.
# 
# In this example, we convert our GeoTIFF file to a NetCDF file. The `-of NetCDF` parameter tells gdal_translate in which format to do the projection. The output file will be named `CMC_glb_TMP.nc`

# In[27]:


get_ipython().run_cell_magic('bash', '', '\ngdal_translate -of NetCDF CMC_glb_TMP.tif CMC_glb_TMP.nc\n')


# Then using `gdalinfo` we can make sure the output NetCDF file is a valid raster file.

# In[29]:


get_ipython().run_cell_magic('bash', '', '\ngdalinfo CMC_glb_TMP.nc\n')


# ## Get the value for a specific point based on a location in longitude/latitude
# 
# Using `gdallocationinfo` command we can get the raw value of a pixel by specifying a location in either longitude/latitude or by specifying a pixel position. 
# 
# In the following example, we use longitude/latitude. The resulting value is in °C.

# In[32]:


get_ipython().run_cell_magic('bash', '', '\ngdallocationinfo -wgs84 CMC_glb_TMP.tif -100 50\n')