<jupyter_start><jupyter_text>------- In this notebook:-------- How to use the Python GDAL extensions to open NetCDF files.- How to reproject datasets between different coordinate systems.- How to combine and merge data from different sources.-------- Required libraries- Himawari Data- Intro to the GDAL library- Digital Elevation model- Data Fusion- Practice these concepts  Canberra, the mountains or the south coast? Where's the best weather?My aunt lives in Narooma, my brother is in Thredbo and I live in Canberra. Last September we interchanged lots of messages in which we were complaining about the weather in our respective places... But, where was the weather the worst? Let's see if we can find out...  0.- To run this notebook on the VDI you'll need to first run:```$ module purge$ module load proj/4.8.0 python/3.5.2 python/3.5.2-matplotlib numpy/1.11.1-py3.5 gdal-python/1.11.1-py3.5 ipython/4.2.0-py3.5$ jupyter notebook```  1.- Intro to Himawari8Himawari 8 is a geostationary satellite that covers the area of the Asia-Pacific and takes a new picture every 10 minutes.The Bureau of Meteorology uses this information operationally for their weather forecasts and the historical data is stored at the NCI at: /g/data2/rr5/satellite/obs/himawari8/FLDK/  1.1.- An index for Himawari8 files<jupyter_code>import datetime
import os.path

def himawari8_path(UTC_datetime, band):
    base = "/g/data2/rr5/satellite/obs/himawari8/FLDK/{year:d}/{month:02d}/{day:02d}/{hour:02d}{minute:02d}/"
    file_name = "{year:d}{month:02d}{day:02d}{hour:02d}{minute:02d}00-P1S-ABOM_BRF_B{band:02d}-PRJ_GEOS141_2000-HIMAWARI8-AHI.nc"

    st = datetime.datetime.strptime(UTC_datetime, "%Y%m%d%H%M")

    path = (base.format(year=st.year, month=st.month, day=st.day, hour=st.hour, minute=st.minute) +
            file_name.format(year=st.year, month=st.month, day=st.day, hour=st.hour, minute=st.minute, band=band))

    if not os.path.isfile(path):
        return None
    
    return path

# For example, to get the path to band 1 taken of the 25th Dec 2015 at 00:00 UTC
himawari8_path("201512250000", 1)<jupyter_output><empty_output><jupyter_text>1.2.- Plot raw data from a Himawari8 file<jupyter_code>%matplotlib inline

from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np

ds = gdal.Open(himawari8_path("201512250000", 1))
band1 = ds.GetRasterBand(1)
raster = band1.ReadAsArray()
no_data = band1.GetNoDataValue()
    
masked = np.ma.array(raster, mask=(raster==no_data))
plt.imshow(masked)<jupyter_output><empty_output><jupyter_text>1.3.- Create a RGB composite using three bands<jupyter_code>def him_rgb(UTC_datetime, rgb_bands, clip):
    red_path = himawari8_path(UTC_datetime, rgb_bands[0])
    ds_r = gdal.Open(red_path)
    red = ds_r.GetRasterBand(1).ReadAsArray()
    red = (red.clip(0, clip) / clip * 255).astype(np.uint8)

    green_path = himawari8_path(UTC_datetime, rgb_bands[1])
    ds_g = gdal.Open(green_path)
    green = ds_g.GetRasterBand(1).ReadAsArray()
    green = (green.clip(0, clip) / clip * 255).astype(np.uint8)

    blue_path = himawari8_path(UTC_datetime, rgb_bands[2])
    ds_b = gdal.Open(blue_path)
    blue = ds_b.GetRasterBand(1).ReadAsArray()
    blue = (blue.clip(0, clip) / clip * 255).astype(np.uint8)

    return np.stack((red, green, blue), axis=2)

rgb = him_rgb("201512250000", [3, 2, 1], .65)
print(rgb.shape)
plt.imshow(rgb)<jupyter_output>(5500, 5500, 3)<jupyter_text>1.4.- Plotting a portion of the image<jupyter_code>plt.imshow(rgb[3200:5000, 1000:3500, :])<jupyter_output><empty_output><jupyter_text>2.- Intro to GDAL: Different datasets cover different areas with different projections and resolutions... How can we combine data coming from different sources?  2.1.- Defining a region in GeoJSON[geojson.io](http://geojson.io)This website allows us to draw polygons and get their gjson representation. For example:'{"type":"Polygon","coordinates":[[[148.5,-36],[148.5,-35],[150.2,-35],[150.2,-36],[148.5,-36]]]}'Represents a box covering the Canberra region. ![caption](./images/notebook_geojson.png)  2.2.- Defining a canvas raster in GDAL: Extent, Size and Projection<jupyter_code>import json

wgs84_wkt ='GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'

def get_raster(geojson_poly, pix_size):
    coords = json.loads(geojson_poly)["coordinates"]
    
    min_lat = min(coords[0][1][1], coords[0][0][1])
    max_lat = max(coords[0][1][1], coords[0][0][1])
    
    min_lon = min(coords[0][2][0], coords[0][0][0])
    max_lon = max(coords[0][2][0], coords[0][0][0])
    
    geotrans = [min_lon, pix_size, 0.0, max_lat, 0.0, -1.0*pix_size]
    size_x = int((max_lon - min_lon) / pix_size)
    size_y = int((max_lat - min_lat) / pix_size)

    mem_drv = gdal.GetDriverByName('MEM')
    ds = mem_drv.Create('', size_x, size_y, 1, gdal.GDT_Float32)
    ds.SetProjection(wgs84_wkt)
    ds.SetGeoTransform(geotrans)

    return ds<jupyter_output><empty_output><jupyter_text>2.3.- Reprojecting between reference systems.<jupyter_code>def reproject(src_file, dst_ds):
    src_ds = gdal.Open(src_file)
    gdal.ReprojectImage(src_ds, dst_ds, src_ds.GetProjection(), dst_ds.GetProjection(), gdal.GRA_NearestNeighbour)

    return dst_ds.GetRasterBand(1).ReadAsArray()<jupyter_output><empty_output><jupyter_text>2.3.- Displaying Himawari 8 for a geoJSON defined area.<jupyter_code>canberra_gjson = '{"type":"Polygon","coordinates":[[[148.5,-36],[148.5,-35],[150.2,-35],[150.2,-36],[148.5,-36]]]}'
him_dest = get_raster(canberra_gjson, .01)
him_proj = reproject(himawari8_path("201512250000", 1), him_dest)

plt.imshow(him_proj, cmap="gray")<jupyter_output><empty_output><jupyter_text>2.4.- A most probably too simplistic way of getting cloud coverage...<jupyter_code>clouds_proj = him_proj > .30
plt.imshow(clouds_proj, cmap="gray")<jupyter_output><empty_output><jupyter_text>2.5.- What is the mean cloud coverage for that image?<jupyter_code>clouds_proj.mean()<jupyter_output><empty_output><jupyter_text>3.- Intro to Digital Elevation Model (DEM) dataDEM models produce regular grid in which each grid point represents a point in a surface and contains the altitude of that location.Geoscience Australia publishes DEM data for Australia with a precission of a seconds of degree (30 meters). This data is stored at the NCI at: /g/data1/rr1/Elevation/  3.1.- Using the above functions to subset these data into the same defined region, resolution and reference system.<jupyter_code>dem_file = "/g/data1/rr1/Elevation/NetCDF/1secSRTM_DEMs_v1.0/DEM-S/Elevation_1secSRTM_DEMs_v1.0_DEM-S_Mosaic_dems1sv1_0.nc"

elv_dest = get_raster(canberra_gjson, .01)
elv_proj = reproject(dem_file, elv_dest)

plt.imshow(elv_proj)<jupyter_output><empty_output><jupyter_text>3.2.- An altitude profile of our selected region:<jupyter_code>plt.hist(elv_proj.ravel(), bins=64, range=(0.0, 2000))
plt.xlabel('Height')
plt.ylabel('Number of Pixels')
plt.show()<jupyter_output><empty_output><jupyter_text>3.3.- Defining our regions of interest:<jupyter_code>elv_canberra = elv_proj[10:50, 70:120]
elv_mountains = elv_proj[50:, :60]
elv_coast = elv_proj[50:, 130:]

print("The average height in Canberra is {0:.2f} m.".format(elv_canberra.mean()))
print("The average height in the mountains is {0:.2f} m.".format(elv_mountains.mean()))
print("The average height in the coast is {0:.2f} m.".format(elv_coast.mean()))<jupyter_output>The average height in Canberra is 776.15 m.
The average height in the mountains is 1245.20 m.
The average height in the coast is 265.97 m.<jupyter_text>4.- Back to out original question: Where is the best weather?  4.1.- Packing cloud coverage in a function<jupyter_code>def give_me_cloud_coverage(utc_datetime, raster_dest):

    him_proj = reproject(himawari8_path(utc_datetime, 1), raster_dest)
    clouds_proj = him_proj > .30
    clouds_canberra = clouds_proj[10:50, 70:120]
    clouds_mountains = clouds_proj[50:, :60]
    clouds_coast = clouds_proj[50:, 130:]

    return int(clouds_canberra.mean()*100), int(clouds_mountains.mean()*100), int(clouds_coast.mean()*100)

cbr, mtn, cst = give_me_cloud_coverage("201609010000", him_dest)

print("Cloud cover 1st Sept 2016. Canberra: {0:d}% Mountains: {1:d}% Coast: {2:d}%".format(cbr, mtn, cst))<jupyter_output>Cloud cover 1st Sept 2016. Canberra: 5% Mountains: 52% Coast: 0%<jupyter_text>4.2.- Collecting cloud data over a period<jupyter_code>def clouds_in_a_period(UTC_start, UTC_end, raster_dest):
    start = datetime.datetime.strptime(UTC_start, "%Y%m%d%H%M")
    end = datetime.datetime.strptime(UTC_end, "%Y%m%d%H%M")
    
    cbr_list = []
    mtn_list = []
    cst_list = []
    while start <= end:
        cbr, mtn, cst = give_me_cloud_coverage(start.strftime("%Y%m%d%H%M"), raster_dest)
        cbr_list.append(cbr)
        mtn_list.append(mtn)
        cst_list.append(cst)
        start = start + datetime.timedelta(hours=24)

    return cbr_list, mtn_list, cst_list

cbr_month, mtn_month, cst_month = clouds_in_a_period("201609010000", "201609300000", him_dest)
days = np.arange(1, 31, 1)
plt.plot(days, cbr_month, 'ro-', days, mtn_month, 'g^-', days, cst_month, 'bs-', label=["a", "b", "c"])
plt.xlabel('Day of the month')
plt.ylabel('Cloud Cover')
lgnd = ['Canberra', 'Mountains', 'Coast']
plt.legend(lgnd)
plt.show()<jupyter_output><empty_output><jupyter_text>4.3.- The veredict<jupyter_code>print("Canberra had an average {0:d}% of clouds during September 2016.".format(int(np.array(cbr_month).mean())))
print("The mountains had and average {0:d}% of clouds during September 2016.".format(int(np.array(mtn_month).mean())))
print("The coast had and average {0:d}% of clouds during September 2016.".format(int(np.array(cst_month).mean())))<jupyter_output>Canberra had an average 57% of clouds during September 2016.
The mountains had and average 66% of clouds during September 2016.
The coast had and average 36% of clouds during September 2016.
{'fraction_non_alphanumeric': 0.13988301774561315, 'fraction_numerical': 0.058094577178546644, 'mean_word_length': 5.258064516129032, 'pattern_counts': {'":': 4, '<': 58, '<?xml version=': 0, '>': 59, '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': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '320240', 'license': 'Apache-2.0', 'n_tokens_mistral': 4171, 'n_tokens_neox': 3763, 'n_words': 862, 'path': 'training/py3/Satellite_Imaging.ipynb', 'repo_name': 'adamsteer/nci-notebooks', 'size': '10087'}