CMAP monthly enhanced mean precipitation data, 1979 - near present

Author: Matt Bertrand

ansible/roles/dataqs/tasks/main.yml

@@ -108,6 +108,11 @@
   command: "{{ app_code_dir }}/venvs/geonode/bin/python {{app_code_dir}}/venvs/geonode/src/dataqs/dataqs/hifld/hifld.py"
   ignore_errors: yes
 
+- include: geoserver_permissions.yml
+
+- name: Create the cmap layer
+  command: "{{ app_code_dir }}/venvs/geonode/bin/python {{app_code_dir}}/venvs/geonode/src/dataqs/dataqs/cmap/cmap.py"
+  ignore_errors: yes
 
 - name: Django updatelayers
   django_manage: command=updatelayers

ansible/roles/dataqs/templates/dataq_settings.py

@@ -35,6 +35,7 @@
     'dataqs.wqp',
     'dataqs.hifld',
     'dataqs.gistemp',
+    'dataqs.cmap',
 )
 
 # CELERY SETTINGS
@@ -127,6 +128,11 @@
         'task': 'dataqs.gistemp.tasks.gistemp_task',
         'schedule': crontab(day_of_month=15, hour=12, minute=0),
         'args': ()
+    },
+    'cmap': {
+        'task': 'dataqs.cmap.tasks.cmap_task',
+        'schedule': crontab(day_of_month=15, hour=12, minute=0),
+        'args': ()
     }
 }
 

dataqs/cmap/__init__.py

@@ -0,0 +1,157 @@
+from __future__ import absolute_import
+import logging
+import os
+import re
+import shutil
+from datetime import date
+from ftplib import FTP
+
+from dateutil.relativedelta import relativedelta
+from dataqs.processor_base import GeoDataMosaicProcessor, GS_DATA_DIR, \
+    GS_TMP_DIR
+from dataqs.helpers import get_band_count, gdal_translate, \
+    nc_convert, style_exists, cdo_fixlng
+
+logger = logging.getLogger("dataqs.processors")
+script_dir = os.path.dirname(os.path.realpath(__file__))
+
+
+class CMAPProcessor(GeoDataMosaicProcessor):
+    """
+    Processor for Land-Ocean Temperature Index, ERSSTv4, 1200km smoothing
+    from the NASA Goddard Institute for Space Studies' Surface Temperature
+    Analysis (GISTEMP).
+    More info at http://data.giss.nasa.gov/gistemp/
+    """
+    prefix = "cmap"
+    base_url = "ftp.cdc.noaa.gov"
+    base_path = "/Datasets/cmap/enh/"
+    base_name = "precip.mon.mean.nc"
+    layer_name = 'cmap'
+    bounds = '-178.75,178.75,-88.75,88.75'
+    title = 'CPC Merged Analysis of Precipitation, 1979/01 - {}'
+    abstract = """The CPC Merged Analysis of Precipitation ("CMAP") is a
+technique which produces pentad and monthly analyses of global precipitation
+ in which observations from raingauges are merged with precipitation estimates
+from several satellite-based algorithms (infrared and microwave). The analyses
+are are on a 2.5 x 2.5 degree latitude/longitude grid and extend back to 1979.
+These data are comparable (but should not be confused with) similarly combined
+analyses by the Global Precipitation Climatology Project which are described in
+Huffman et al (1997).
+
+It is important to note that the input data sources to make these analyses are
+not constant throughout the period of record. For example, SSM/I (passive
+microwave - scattering and emission) data became available in July of 1987;
+prior to that the only microwave-derived estimates available are from the MSU
+algorithm (Spencer 1993) which is emission-based thus precipitation estimates
+are avaialble only over oceanic areas. Furthermore, high temporal resolution IR
+data from geostationary satellites (every 3-hr) became available during 1986;
+prior to that, estimates from the OPI technique (Xie and Arkin 1997) are used
+based on OLR from polar orbiting satellites.
+
+The merging technique is thoroughly described in Xie and Arkin (1997). Briefly,
+the methodology is a two-step process. First, the random error is reduced by
+linearly combining the satellite estimates using the maximum likelihood method,
+in which case the linear combination coefficients are inversely proportional to
+the square of the local random error of the individual data sources. Over global
+ land areas the random error is defined for each time period and grid location
+by comparing the data source with the raingauge analysis over the surrounding
+area. Over oceans, the random error is defined by comparing the data sources
+with the raingauge observations over the Pacific atolls. Bias is reduced when
+the data sources are blended in the second step using the blending technique of
+Reynolds (1988). Here the data output from step 1 is used to define the "shape"
+of the precipitation field and the rain gauge data are used to constrain the
+amplitude.
+
+Monthly and pentad CMAP estimates back to the 1979 are available from CPC ftp
+server.
+
+References:
+
+Huffman, G. J. and co-authors, 1997: The Global Precipitation Climatology
+Project (GPCP) combined data set. Bull. Amer. Meteor. Soc., 78, 5-20.
+
+Reynolds, R. W., 1988: A real-time global sea surface temperature analysis. J.
+Climate, 1, 75-86.
+
+Spencer, R. W., 1993: Global oceanic precipitation from the MSU during 1979-91
+and comparisons to other climatologies. J. Climate, 6, 1301-1326.
+
+Xie P., and P. A. Arkin, 1996: Global precipitation: a 17-year monthly analysis
+based on gauge observations, satellite estimates, and numerical model outputs.
+Bull. Amer. Meteor. Soc., 78, 2539-2558.
+"""
+
+    def download(self, url, tmp_dir=GS_TMP_DIR, filename=None):
+        if not filename:
+            filename = url.rsplit('/')[-1]
+        ftp = FTP(self.base_url)
+        ftp.login('anonymous', 'anonymous')
+        ftp.cwd(self.base_path)
+        with open(os.path.join(self.tmp_dir, filename),
+                  'wb') as outfile:
+            ftp.retrbinary('RETR %s' % self.base_name, outfile.write)
+        return filename
+
+    def convert(self, nc_file):
+        nc_transform = nc_convert(nc_file)
+        cdo_transform = cdo_fixlng(nc_transform, bounds=self.bounds)
+        return cdo_transform
+
+    def extract_band(self, tif, band, outname):
+        outfile = os.path.join(self.tmp_dir, outname)
+        gdal_translate(tif, outfile, bands=[band],
+                       projection='EPSG:4326',
+                       options=['TILED=YES', 'COMPRESS=LZW'])
+        return outfile
+
+    def get_date(self, months):
+        start_month = date(1979, 1, 1)
+        return start_month + relativedelta(months=months - 1)
+
+    def get_title(self, months):
+        end_month = self.get_date(months)
+        return self.title.format(end_month.strftime('%Y/%m'))
+
+    def run(self):
+        """
+        Retrieve and process the latest NetCDF file.
+        """
+        ncfile = self.download(
+            self.base_url, filename='{}.nc'.format(self.layer_name))
+        cdf_file = self.convert(os.path.join(self.tmp_dir, ncfile))
+        bands = get_band_count(cdf_file)
+        img_list = self.get_mosaic_filenames(self.layer_name)
+        for band in range(1, bands + 1):
+            band_date = re.sub('[\-\.]+', '', self.get_date(band).isoformat())
+            img_name = '{}_{}T000000000Z.tif'.format(self.layer_name, band_date)
+            if img_name not in img_list:
+                band_tif = self.extract_band(cdf_file, band, img_name)
+                dst_file = self.data_dir.format(gsd=GS_DATA_DIR,
+                                                ws=self.workspace,
+                                                layer=self.layer_name,
+                                                file=img_name)
+                dst_dir = os.path.dirname(dst_file)
+                if not os.path.exists(dst_dir):
+                    os.makedirs(dst_dir)
+                if dst_file.endswith('.tif'):
+                    shutil.move(os.path.join(self.tmp_dir, band_tif), dst_file)
+                    self.post_geoserver(dst_file, self.layer_name)
+
+        if not style_exists(self.layer_name):
+            with open(os.path.join(script_dir,
+                                   'resources/cmap.sld')) as sld:
+                self.set_default_style(self.layer_name, self.layer_name,
+                                       sld.read().format(latest_band=bands))
+        self.update_geonode(self.layer_name, title=self.get_title(bands),
+                            description=self.abstract,
+                            store=self.layer_name,
+                            bounds=('-178.75', '178.75', '-88.75', '88.75',
+                                    'EPSG:4326'))
+        self.truncate_gs_cache(self.layer_name)
+        self.cleanup()
+
+
+if __name__ == '__main__':
+    processor = CMAPProcessor()
+    processor.run()

dataqs/cmap/cmap.py

@@ -0,0 +1,32 @@
+<?xml version="1.0" encoding="UTF-8"?><sld:StyledLayerDescriptor xmlns="http://www.opengis.net/sld" xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" version="1.0.0">
+  <sld:NamedLayer>
+    <sld:Name>cmap</sld:Name>
+    <sld:UserStyle>
+      <sld:Name>cmap</sld:Name>
+      <sld:Title>cmap</sld:Title>
+      <sld:FeatureTypeStyle>
+        <sld:Name>name</sld:Name>
+        <sld:Rule>
+          <sld:RasterSymbolizer>
+        <Opacity>1.0</Opacity>
+        <ChannelSelection>
+                <GrayChannel>
+                        <SourceChannelName>{latest_band}</SourceChannelName>
+                </GrayChannel>
+        </ChannelSelection>
+        <ColorMap extended="true">
+            <sld:ColorMapEntry color="#2b83ba" label="0 mm" opacity="1.0" quantity="0"/>
+            <sld:ColorMapEntry color="#74b6ad" label="10 mm" opacity="1.0" quantity="10"/>
+            <sld:ColorMapEntry color="#b7e1a7" label="20 mm" opacity="1.0" quantity="20"/>
+            <sld:ColorMapEntry color="#e7f5b7" label="30 mm" opacity="1.0" quantity="30"/>
+            <sld:ColorMapEntry color="#fee7a4" label="40 mm" opacity="1.0" quantity="40"/>
+            <sld:ColorMapEntry color="#fdb96e" label="50 mm" opacity="1.0" quantity="50"/>
+            <sld:ColorMapEntry color="#ec6e43" label="60 mm" opacity="1.0" quantity="60"/>
+            <sld:ColorMapEntry color="#d7191c" label="70 mm" opacity="1.0" quantity="70"/>
+        </ColorMap>
+          </sld:RasterSymbolizer>
+        </sld:Rule>
+      </sld:FeatureTypeStyle>
+    </sld:UserStyle>
+  </sld:NamedLayer>
+</sld:StyledLayerDescriptor>

dataqs/cmap/resources/cmap.nc

@@ -0,0 +1,10 @@
+from __future__ import absolute_import
+
+from celery import shared_task
+from dataqs.cmap.cmap import CMAPProcessor
+
+
+@shared_task
+def cmap_task():
+    processor = CMAPProcessor()
+    processor.run()

dataqs/cmap/resources/cmap.sld

@@ -0,0 +1,70 @@
+import glob
+import os
+import unittest
+from datetime import date
+from django.test import TestCase
+from dataqs.cmap.cmap import CMAPProcessor
+
+script_dir = os.path.dirname(os.path.realpath(__file__))
+
+
+class CMAPPTest(TestCase):
+    """
+    Tests the dataqs.cmap module.  Since each processor is highly
+    dependent on a running GeoNode instance for most functions, only
+    independent functions are tested here.
+    """
+
+    def setUp(self):
+        self.processor = CMAPProcessor()
+        self.dl_file = os.path.join(script_dir, 'resources/cmap.nc')
+
+    def tearDown(self):
+        self.processor.cleanup()
+
+    def test_convert(self):
+        """
+        Verify that a NetCDF file is converted via netcdf and cdo apps.
+        Skip if these apps are not installed.
+        """
+        try:
+            converted_nc = self.processor.convert(self.dl_file)
+            self.assertTrue(os.path.exists(converted_nc))
+        except OSError:
+            # cdo and/or netcdf not installed
+            raise unittest.SkipTest()
+
+    def test_extract_band(self):
+        """
+        Verify that a GeoTIFF file is created.
+        """
+        dl_tif = self.processor.extract_band(
+            self.dl_file, 1, os.path.join(self.processor.tmp_dir, 'cmap.tif'))
+        self.assertTrue(os.path.exists(dl_tif))
+
+    def test_get_title(self):
+        """
+        Verify that the correct title is returned
+        """
+        title = self.processor.get_title(451)
+        self.assertEquals(
+            'CPC Merged Analysis of Precipitation, 1979/01 - 2016/07', title)
+
+    def test_get_date(self):
+        """
+        Verify that the correct date is returned
+        """
+        band_date = self.processor.get_date(451)
+        self.assertEquals(band_date, date(2016, 7, 1))
+
+    def test_cleanup(self):
+        """
+        Make sure temporary files are deleted.
+        """
+        self.processor.extract_band(
+            self.dl_file, 1, os.path.join(self.processor.tmp_dir, 'cmap.tif'))
+        self.assertNotEqual([], glob.glob(os.path.join(
+            self.processor.tmp_dir, self.processor.prefix + '*')))
+        self.processor.cleanup()
+        self.assertEquals([], glob.glob(os.path.join(
+            self.processor.tmp_dir, self.processor.prefix + '*')))

dataqs/cmap/tasks.py

@@ -124,7 +124,7 @@ def nc_convert(filename):
     """
     output_file = re.sub('\.nc$', '.classic.nc', filename)
     subprocess.check_call(
-        ["/usr/local/bin/nccopy", "-k", "classic", filename, output_file])
+        ["nccopy", "-k", "classic", filename, output_file])
     return output_file
 
 
@@ -143,15 +143,16 @@ def cdo_invert(filename):
     return output_file
 
 
-def cdo_fixlng(filename):
+def cdo_fixlng(filename, bounds="-180,180,-90,90"):
     """
     Change the longitude coordinates of a NetCDF file from 0->360 to -180->180
     :param filename: Full path * name of NetCDF image to process
+    :param bounds: String representation of bounds (minX,maxX,minY,maxY)
     :return: output filename
     """
 
     output_file = "{}.lng.nc".format(os.path.splitext(filename)[0])
-    subprocess.check_call(["cdo", "sellonlatbox,-180,180,-90,90", "{}".format(
+    subprocess.check_call(["cdo", "sellonlatbox,{}".format(bounds), "{}".format(
         filename), output_file])
     return output_file