Use Radar Data from Obs Via Amazon S3 bucket

Author: wxtim

cylc/flow/etc/tutorial/cylc-forecasting-workflow/bin/forecast

@@ -150,11 +150,12 @@ def push_rainfall(rainfall, wind_data, step, resolution, spline_level):
                                 dim_x, dim_y, resolution, resolution,
                                 spline_level)
 
+    domain = util.parse_domain(os.environ['DOMAIN'])
+
     while True:
         out_of_bounds = []
         for itt in range(len(x_values)):
             try:
-                domain = util.parse_domain(os.environ['DOMAIN'])
                 lng = domain['lng1'] + x_values[itt]
                 lat = domain['lat1'] + y_values[itt]
 

cylc/flow/etc/tutorial/cylc-forecasting-workflow/bin/get-observations

@@ -73,19 +73,6 @@ class NoDataException(Exception):
     ...
 
 
-def get_datapoint_data(site_id, time, api_key):
-    """Get weather data from the DataPoint service."""
-    # The URL required to get the data.
-    time = datetime.strptime(time, '%Y%m%dT%H%MZ').strftime('%Y-%m-%dT%H:%MZ')
-    url = BASE_URL.format(time=time, site_id=site_id, api_key=api_key)
-    req = requests.get(url)
-    if req.status_code != 200:
-        raise Exception(f'{url} returned exit code {req.status_code}')
-    # Get the data and parse it as JSON.
-    print('Opening URL: %s' % url)
-    return req.json()['SiteRep']['DV']['Location']
-
-
 def get_archived_data(site_id, time):
     """Return archived weather data from the workflow directory."""
     print(
@@ -202,13 +189,8 @@ def main(site_id, api_key=None):
             site_id, dist = get_nearby_site(site_id, badsites)
 
     if obs is None:
-        if api_key:
-            print('Attempting to get weather data from DataPoint...')
-            data = get_datapoint_data(site_id, cycle_point, api_key)
-        else:
-            print('No API key provided, falling back to archived data')
-            data = get_archived_data(site_id, cycle_point)
-
+        print('No API key provided, falling back to archived data')
+        data = get_archived_data(site_id, cycle_point)
         obs = extract_observations(data)
 
     # Write observations.

cylc/flow/etc/tutorial/cylc-forecasting-workflow/bin/get-rainfall

@@ -38,20 +38,20 @@ import math
 import os
 import shutil
 
+import h5py
+from pathlib import Path
+import urllib
 import requests
 
-try:
-    from PIL import Image
-except ModuleNotFoundError:
-    # not all PIL installations are created equal
-    # sometimes we must import like this
-    import Image
 from mercator import get_offset, get_scale, pos_to_coord
 import util
 
-
-URL = ('http://datapoint.metoffice.gov.uk/public/data/layer/wxobs/'
-       'RADAR_UK_Composite_Highres/png?TIME={time}&key={api_key}')
+S3URL = (
+    'https://met-office-radar-obs-data.s3-eu-west-2.amazonaws.com/radar/'
+    '{Y}/{m}/{d}/{YYYYmmddHHMM}_ODIM_ng_radar_rainrate_composite_1km_UK.h5'
+)
+DEBUG=os.environ['CYLC_DEBUG']=='true'
+CYLC_TASK_LOG_DIR=os.environ['CYLC_TASK_LOG_DIR']
 
 
 class Rainfall(object):
@@ -62,16 +62,6 @@ class Rainfall(object):
         resolution (float): The length of each grid cell in degrees.
 
     """
-    VALUE_MAP = {
-        (0, 0, 254, 255): 1,
-        (50, 101, 254, 255): 2,
-        (127, 127, 0, 255): 3,
-        (254, 203, 0, 255): 4,
-        (254, 152, 0, 255): 5,
-        (254, 0, 0, 255): 6,
-        (254, 0, 254, 255): 7
-    }
-
     def __init__(self, domain, resolution):
         self.resolution = resolution
         self.domain = domain
@@ -98,10 +88,23 @@ class Rainfall(object):
         """
         itt_x, itt_y = util.get_grid_coordinates(lng, lat, self.domain,
                                                  self.resolution)
-        try:
-            self.data[itt_y][itt_x].append(self.VALUE_MAP[value])
-        except KeyError:
-            pass
+                    
+        self.data[itt_y][itt_x].append(self.value_map(value))
+
+    @staticmethod
+    def value_map(v):
+        """Convert rainfall rate values into colour space values.
+
+        Checks if rainfall value above each threshold in turn.
+
+        TODO:
+            - Unit test this
+        """
+        thresholds = {32, 16, 8, 4, 2, 1, .5, .2}
+        for i, threshold in enumerate(sorted(thresholds, reverse=True)):
+            if v > threshold:
+                return 8 - i
+        return 0
 
     def compute_bins(self):
         """Return this dataset as a 2D matrix."""
@@ -114,25 +117,6 @@ class Rainfall(object):
         return self.data
 
 
-def get_datapoint_radar_image(filename, time, api_key):
-    """Retrieve a png image of rainfall from the DataPoint service.
-
-    Args:
-        filename (str): The path to write the image file to.
-        time (str): The datetime of the image to retrieve in ISO8601 format.
-        api_key (str): Datapoint API key.
-
-    """
-    time = datetime.strptime(time, '%Y%m%dT%H%MZ').strftime(
-        '%Y-%m-%dT%H:%M:%SZ')
-    url = URL.format(time=time, api_key=api_key)
-    req = requests.get(url)
-    if req.status_code != 200:
-        raise Exception(f'{url} returned exit code {req.status_code}')
-    with open(filename, 'bw') as png_file:
-        png_file.write(req.content)
-
-
 def get_archived_radar_image(filename, time):
     """Retrieve a png image from the archived data in the workflow directory.
 
@@ -147,8 +131,23 @@ def get_archived_radar_image(filename, time):
         filename)
 
 
+def get_amazon_radar_data(filename, time):
+    time = datetime.strptime(time, '%Y%m%dT%H%MZ')
+    url = S3URL.format(
+        Y=time.strftime('%Y'),
+        m=time.strftime('%m'),
+        d=time.strftime('%d'),
+        YYYYmmddHHMM=time.strftime('%Y%m%d%H%M')
+    )
+    print(f'[INFO] Getting data from {url=}')
+    data = urllib.request.urlopen(url).read()
+    with open(filename, 'wb') as fh:
+        fh.write(data)
+
+
+# def process_rainfall_data(filename, resolution, domain):
 def process_rainfall_data(filename, resolution, domain):
-    """Generate a 2D matrix of data from the rainfall data in the image.
+    """get_amazon_radar_dataGenerate a 2D matrix of data from the rainfall data in the image.
 
     Args:
         filename (str): Path to the png image to process.
@@ -160,37 +159,64 @@ def process_rainfall_data(filename, resolution, domain):
         list - A 2D matrix of rainfall data.
 
     """
+    print(f'{util.INFO}Analysing data from {filename}')
+    data = h5py.File(filename)['dataset1']['data1']['data']
     rainfall = Rainfall(domain, resolution)
 
-    image = Image.open(filename)
-    scale = get_scale(domain, image.width)
+    # image = Image.open(filename)
+    height, width = data.shape
+
+    # TODO - work out why this isn't working
+    scale = get_scale(domain, width)
     offset = get_offset(domain, scale)
+    scale = (1672.2334443981335, 3344.466888796267)
+    offset = (-1100.8461538461539, 1400.6953225710452)
+
+    if DEBUG:
+        from matplotlib import pyplot as plt
+        Path('CYLC_TASK_LOG_DIR').mkdir(parents=True, exist_ok=True)
+        plt.imshow(data)
+        plt.savefig(f'{os.environ['CYLC_TASK_LOG_DIR']}/raw.data.png')
 
-    for itt_x in range(image.width):
-        for itt_y in range(image.height):
+    for itt_x in range(width):
+        for itt_y in range(height):
             lng, lat = pos_to_coord(
                 itt_x,
                 itt_y * (2. / 3.),  # Counter aspect ratio.
-                offset, scale)
-            rainfall.add(lng, lat, image.getpixel((itt_x, itt_y)))
+                offset,
+                scale
+            )
+            val = float(data[height - 1 - itt_y][itt_x])
+            val = 0 if val == -1 else val
+            rainfall.add(lng, lat, val)
+    data = rainfall.compute_bins()
 
-    return rainfall.compute_bins()
+    if DEBUG:
+        plt.imshow(data)
+        plt.savefig(f'{CYLC_TASK_LOG_DIR}/processed.data.png')
 
+    return data
 
 def main():
     time = os.environ['CYLC_TASK_CYCLE_POINT']
     resolution = float(os.environ['RESOLUTION'])
     domain = util.parse_domain(os.environ['DOMAIN'])
-    api_key = os.environ.get('API_KEY')
 
-    if api_key:
-        print('Attempting to get weather data from the DataPoint service.')
-        get_datapoint_radar_image('rainfall-radar.png', time, api_key)
+    # Acts as a switch - if a file-name is given, use that file-name
+    # TODO - keep the nice S3 formatting with implied metadata
+    canned_data = os.environ.get('CANNED_DATA', 'fetch')
+
+    if canned_data == 'fetch':
+        print(f'{util.INFO}Attempting to get rainfall data from S3 bucket')
+        get_amazon_radar_data('radardata.h5', time)
+        canned_data = 'radardata.h5'
     else:
-        print('No API key provided, falling back to archived data')
-        get_archived_radar_image('rainfall-radar.png', time)
+        print(
+            f'{util.INFO}Canned data provided: {canned_data}')
+        # get_archived_radar_data(canned_data, time)
+
+    data = process_rainfall_data(canned_data, resolution, domain)
 
-    data = process_rainfall_data('rainfall-radar.png', resolution, domain)
     util.write_csv('rainfall.csv', data)
 
 

cylc/flow/etc/tutorial/cylc-forecasting-workflow/flow.cylc

@@ -39,15 +39,24 @@
         [[[environment]]]
             # The dimensions of each grid cell in degrees.
             RESOLUTION = 0.2
+
             # The area to generate forecasts for (lng1, lat1, lng2, lat2)
-            DOMAIN = -12,48,5,61  # Do not change!
+            # Amazon S3 domain (my interpolation)
+            DOMAIN = -14,44,12,63
+
+        [[[meta]]]
+            metadata source = Met Office Observations R&D
+            ellipsiod = "Airy 1830"
+            projection = "Transverse Mercator"
+            origin = "49,2"
+            SW = 43.8341,-12.0118
+            NW = 62.9207,-17.9739
+            NE = 62.6515,16.1289
+            SE = 437009,941578
+            old domain = -12,48,5,61
 
     [[get_observations<station>]]
         script = get-observations
-        [[[environment]]]
-            # The key required to get weather data from the DataPoint service.
-            # To use archived data comment this line out.
-            API_KEY = DATAPOINT_API_KEY
 
     [[get_observations<station=aldergrove>]]
         [[[environment]]]
@@ -70,7 +79,13 @@
         [[[environment]]]
             # The key required to get weather data from the DataPoint service.
             # To use archived data comment this line out.
-            API_KEY = DATAPOINT_API_KEY
+            API_KEY = amazon
+            MAP_FILE = "${CYLC_TASK_LOG_ROOT}-map.html"
+
+            # Create the html map file in the task's log directory.
+            MAP_FILE = "${CYLC_TASK_LOG_ROOT}-map.html"
+            # The path to the template file used to generate the html map.
+            MAP_TEMPLATE = "$CYLC_WORKFLOW_RUN_DIR/lib/template/map.html"
 
     [[forecast]]
         script = forecast 60 5  # Generate 5 forecasts at 60 minute intervals.
@@ -90,5 +105,6 @@
     [[post_process<site>]]
         # Generate a forecast for the location <site> 60 minutes in the future.
         script = post-process $CYLC_TASK_PARAM_site 60
+        post-script = sleep 365000
 
 {% include 'etc/python-job.settings' %}

cylc/flow/etc/tutorial/cylc-forecasting-workflow/lib/python/util.py

@@ -28,6 +28,7 @@
 import jinja2
 
 R_0 = 6371.  # Radius of the Earth (km).
+INFO = '[INFO] '
 
 
 def frange(start, stop, step):