Initial

README.md

@@ -0,0 +1,52 @@
+![WMO stations](img/stations.jpg)
+
+This is a commandline tool to plot global temperature data using data from the Global Historical Climatology Network.
+
+See https://www.ncei.noaa.gov/pub/data/ghcn/v4
+
+Installation
+============
+
+To compile from source first install the prerequisites:
+
+On a Debian based system:
+
+``` bash
+sudo apt-get install build-essential gnuplot wget marble python3
+```
+
+On Arch/Parabola:
+
+``` bash
+sudo pacman -S gnuplot wget marble python
+```
+
+Obtaining the data
+==================
+
+Before you begin you'll need to download the latest version of the GHCN version 4 data.  tempgraph needs three files: the country codes, the weather stations and the temperature data itself.  These files can be obtained here:
+
+```bash
+mkdir data
+cd data
+wget https://www.ncei.noaa.gov/pub/data/ghcn/v4/ghcnm.tavg.latest.qcf.tar.gz
+wget https://www.ncei.noaa.gov/pub/data/ghcn/v4/ghcnm-countries.txt
+tar -xzvf ghcnm.tavg.latest.qcf.tar.gz
+cd ..
+```
+
+The compressed archive contains two files, one which is the temperature data (.dat) and the other which contains details of the weather stations (.inv). I typically rename these to:
+
+```bash
+cp data/ghcnm.v4*/*.dat data/v4.mean
+cp data/ghcnm.v4*/*.inv data/wmo.txt
+cp data/ghcnm-countries.txt data/v4.country.codes
+```
+
+Usage
+=====
+
+``` bash
+python3 tempgraph2.py
+ls *.kml *.jpg
+```

anomaly.py

@@ -0,0 +1,137 @@
+__filename__ = "anomaly.py"
+__author__ = "Bob Mottram"
+__license__ = "GPL3+"
+__version__ = "2.0.0"
+__maintainer__ = "Bob Mottram"
+__email__ = "[email protected]"
+__status__ = "Production"
+__module_group__ = "Commandline Interface"
+
+import os
+
+
+def _stationsAnomaly(year: int, stationsData: {},
+                     baseline: [], stationIds: set) -> float:
+    """Returns anomaly for the given stations for the given year
+    """
+    anomaly = 0.0
+    hits = 0
+
+    for id in stationIds:
+        if id not in stationsData:
+            continue
+        if year not in stationsData[id]:
+            continue
+        monthData = stationsData[id][year]['month']
+        for monthIndex in range(12):
+            if not baseline[monthIndex]:
+                continue
+            if monthData[monthIndex]['av'] > -80:
+                if monthData[monthIndex]['av'] < 80:
+                    anomaly += \
+                        monthData[monthIndex]['av'] - baseline[monthIndex]
+                    hits += 1
+    if hits > 0:
+        return anomaly / float(hits)
+    return None
+
+
+def updateGridAnomalies(grid: [], stationsData: {},
+                        startYear: int, endYear: int) -> int:
+    """Calculates anomalies for each grid cell within a range of years
+    """
+    ctr = 0
+    yearCtr = 0
+    for gridCell in grid:
+        gridCell['anomalies'] = {}
+        if not gridCell['stationIds']:
+            continue
+        baseline = gridCell['baseline']
+        stationIds = gridCell['stationIds']
+        for year in range(startYear, endYear + 1, 1):
+            gridCell['anomalies'][year] = \
+                _stationsAnomaly(year, stationsData, baseline, stationIds)
+            if gridCell['anomalies'][year] is not None:
+                yearCtr += 1
+            ctr += 1
+    if yearCtr > 0:
+        return int(yearCtr * 100 / float(ctr))
+    return 0
+
+
+def getGlobalAnomalies(grid: [],
+                       startYear: int, endYear: int) -> {}:
+    """Returns global anomalies in the given year range
+    """
+    result = {}
+    for year in range(startYear, endYear + 1, 1):
+        yearAnomaly = 0.0
+        ctr = 0
+        for gridCell in grid:
+            if gridCell['anomalies'].get(year):
+                yearAnomaly += gridCell['anomalies'][year]
+                ctr += 1
+        if ctr > 0:
+            result[year] = yearAnomaly / float(ctr)
+        else:
+            result[year] = None
+    return result
+
+
+def plotGlobalAnomalies(grid: [],
+                        startYear: int, endYear: int) -> None:
+    """Plot anomalies graph
+    """
+    anomalies = getGlobalAnomalies(grid, startYear, endYear)
+    series = []
+    minimumTemp = 99999999
+    maximumTemp = -99999999
+    for year in range(startYear, endYear + 1, 1):
+        series.append(anomalies[year])
+        if anomalies[year] > maximumTemp:
+            maximumTemp = anomalies[year]
+        if anomalies[year] < minimumTemp:
+            minimumTemp = anomalies[year]
+
+    title = \
+        "Global Temperature Anomalies " + \
+        str(startYear) + ' - ' + str(endYear)
+    subtitle = "Source https://www.ncei.noaa.gov/pub/data/ghcn/v4"
+    Xlabel = 'Year'
+    Ylabel = 'Average Temperature Anomaly (Celcius)'
+    indent = 0.34
+    vpos = 0.94
+    imageWidth = 1000
+    imageHeight = 1000
+    plotName = 'global_anomalies'
+    imageFormat = 'jpg'
+    imageFormat2 = 'jpeg'
+    filename = plotName + '.' + imageFormat
+    scriptFilename = plotName + '.gnuplot'
+    dataFilename = plotName + '.data'
+    with open(dataFilename, 'w+') as fp:
+        year = startYear
+        for temperatureAnomaly in series:
+            fp.write(str(year) + "    " + str(temperatureAnomaly) + '\n')
+            year += 1
+    script = \
+        "reset\n" + \
+        "set title \"" + title + "\"\n" + \
+        "set label \"" + subtitle + "\" at screen " + \
+        str(indent) + ", screen " + str(vpos) + "\n" + \
+        "set yrange [" + str(minimumTemp) + ":" + \
+        str(maximumTemp) + "]\n" + \
+        "set xrange [" + str(startYear) + ":" + str(endYear) + "]\n" + \
+        "set lmargin 9\n" + \
+        "set rmargin 2\n" + \
+        "set xlabel \"" + Xlabel + "\"\n" + \
+        "set ylabel \"" + Ylabel + "\"\n" + \
+        "set grid\n" + \
+        "set key right bottom\n" + \
+        "set terminal " + imageFormat2 + \
+        " size " + str(imageWidth) + "," + str(imageHeight) + "\n" + \
+        "set output \"" + filename + "\"\n" + \
+        "plot \"" + dataFilename + "\" using 1:2 notitle with lines\n"
+    with open(scriptFilename, 'w+') as fp:
+        fp.write(script)
+    os.system('gnuplot ' + scriptFilename)

baseline.py

@@ -0,0 +1,102 @@
+__filename__ = "baseline.py"
+__author__ = "Bob Mottram"
+__license__ = "GPL3+"
+__version__ = "2.0.0"
+__maintainer__ = "Bob Mottram"
+__email__ = "[email protected]"
+__status__ = "Production"
+__module_group__ = "Commandline Interface"
+
+
+def _getBaselineForYear(id: str, stationsData: {}, year: int) -> []:
+    """Returns the baseline for a given station id
+    """
+    baseline = [None] * 12
+    if year not in stationsData[id]:
+        return baseline
+    monthData = stationsData[id][year]['month']
+    for monthIndex in range(12):
+        monthData2 = monthData[monthIndex]
+        if monthData2['qcflag'] == 'M':
+            # Flagged as error
+            continue
+        if monthData2['av'] > -80:
+            if monthData2['av'] < 80:
+                baseline[monthIndex] = monthData2['av']
+    return baseline
+
+
+def _getBaseline(id: str, stationsData: {},
+                 startYear: int, endYear: int) -> []:
+    """Returns the baseline for a given station id
+    """
+    if id not in stationsData:
+        return [None] * 12
+
+    if startYear == endYear:
+        return _getBaselineForYear(id, stationsData, startYear)
+    else:
+        baseline = [0.0] * 12
+        hits = [0] * 12
+        for year in range(startYear, endYear, 1):
+            if year not in stationsData[id]:
+                continue
+            monthData = stationsData[id][year]['month']
+            for monthIndex in range(12):
+                monthData2 = monthData[monthIndex]
+                if monthData2['qcflag'] == 'M':
+                    # Flagged as error
+                    continue
+                if monthData2['av'] > -80:
+                    if monthData2['av'] < 80:
+                        baseline[monthIndex] += monthData2['av']
+                        hits[monthIndex] += 1
+
+        for monthIndex in range(12):
+            if hits[monthIndex] > 0:
+                baseline[monthIndex] /= float(hits[monthIndex])
+            else:
+                baseline[monthIndex] = None
+    return baseline
+
+
+def _baselineForStations(stationsData: {},
+                         startYear: int, endYear: int,
+                         stationIds: set) -> []:
+    """Returns a baseline for the given range of years
+    and the given station ids
+    """
+    if not stationIds:
+        return [None] * 12, False
+
+    baseline = [0.0] * 12
+    hits = [0] * 12
+    for id in stationIds:
+        stationBaseline = _getBaseline(id, stationsData, startYear, endYear)
+        for monthIndex in range(12):
+            if stationBaseline[monthIndex] is not None:
+                baseline[monthIndex] += stationBaseline[monthIndex]
+                hits[monthIndex] += 1
+
+    for monthIndex in range(12):
+        if hits[monthIndex] > 0:
+            baseline[monthIndex] /= float(hits[monthIndex])
+        else:
+            baseline[monthIndex] = None
+
+    return baseline, True
+
+
+def updateGridBaselines(grid: [], stationsData: {},
+                        startYear: int, endYear: int) -> int:
+    """Calculates reference baselines for each grid cell
+    """
+    ctr = 0
+    for gridCell in grid:
+        if gridCell['stationIds']:
+            gridCell['baseline'], hasData = \
+                _baselineForStations(stationsData, startYear, endYear,
+                                     gridCell['stationIds'])
+            if hasData:
+                ctr += 1
+    return ctr