diff --git a/.gitignore b/.gitignore index 4a3a461..b55a906 100644 --- a/.gitignore +++ b/.gitignore @@ -36,6 +36,9 @@ grabAQfiles.sh Python/__pycache__/ Python/*.pyc +Python/*/__pycache__/ +Python/*/*.pyc +Python/*/*av*.txt CALPUFF_EXE/*.exe CALPUFF_EXE/libraries/lib* diff --git a/Python/analysis_tools/conc.nc b/Python/analysis_tools/conc.nc new file mode 100644 index 0000000..377e0e9 Binary files /dev/null and b/Python/analysis_tools/conc.nc differ diff --git a/Python/analysis_tools/images/Avconc_stamp.png b/Python/analysis_tools/images/Avconc_stamp.png new file mode 100644 index 0000000..ea82dbe Binary files /dev/null and b/Python/analysis_tools/images/Avconc_stamp.png differ diff --git a/Python/analysis_tools/images/P_stamp.png b/Python/analysis_tools/images/P_stamp.png new file mode 100644 index 0000000..126f5c6 Binary files /dev/null and b/Python/analysis_tools/images/P_stamp.png differ diff --git a/Python/analysis_tools/images/test.png b/Python/analysis_tools/images/test.png new file mode 100644 index 0000000..72921fb Binary files /dev/null and b/Python/analysis_tools/images/test.png differ diff --git a/Python/analysis_tools/maptoolkit.py b/Python/analysis_tools/maptoolkit.py new file mode 120000 index 0000000..4c909b0 --- /dev/null +++ b/Python/analysis_tools/maptoolkit.py @@ -0,0 +1 @@ +../maptoolkit.py \ No newline at end of file diff --git a/Python/analysis_tools/monthly.py b/Python/analysis_tools/monthly.py new file mode 100644 index 0000000..745b3a6 --- /dev/null +++ b/Python/analysis_tools/monthly.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""pmaps +.. module:: genmaps + :platform: Unix + :synopis: +.. moduleauther: CEMAC (UoL) +.. description: This module was developed by CEMAC as part of the UNRESP + Project. This Script plots CALPUFF concrec data on a map. + :copyright: © 2019 University of Leeds. + :license: BSD-2 Clause. +Example: + To use:: + ./genemaps.py + - Date string, format YYYYMMDD, of the current CALPUFF run. + Used to locate directory containing the SO2 output files + (with assumed naming convention 'concrec0100**.dat', + where '**' goes from '01' through to '48') +.. CEMAC_UNRESPForcastingSystem: + https://github.com/cemac/UNRESPForcastingSystem +""" +from __future__ import print_function +import argparse +from dateutil.parser import parse +import datetime +from maptoolkit import * + + +def conc_array_plain(ny, nx, filePaths): + # Read in concentration data: + f = open(filePaths, 'r') + lines = f.read().splitlines() + f.close + # Process concentration data into desired format: + conc = np.array([float(X) for X in lines]) * 100**3 # ug/cm^3 -> ug/m^3 + concAry = np.reshape(conc, (ny, nx)) # Reshape data onto latlon grid + return concAry, conc + + +# READ IN COMMAND LINE ARGUMENTS +dstring = ("Date string, format YYYYMMDD, START date and End Date") +parser = argparse.ArgumentParser(description=dstring) +parser.add_argument("--start", help=dstring, type=str) +parser.add_argument("--end", help=dstring, type=str) +args = parser.parse_args() + +# Required files and info +start = datetime.datetime.strptime(args.start, "%Y%m%d") +end = datetime.datetime.strptime(args.end, "%Y%m%d") +datahome = "/nfs/see-fs-01_users/earhbu/UNRESP_SPACE/UNRESPForecastingSystem/" +nConcFiles = 48 +binLims = [10, 350, 600, 2600, 9000, 14000] # SO2 bin limits +xyFile = "../../data/xy_masaya.dat" +glat, glon, latMin, latMax, lonMin, lonMax, ny, nx = genxy(xyFile) +datelist = [start + datetime.timedelta(days=x) for x in range(0, (end-start).days)] +nhours = 0 +concsum = np.zeros(18000) +for date in datelist: + month = str(date.strftime("%Y%m")) + year = str(date.strftime("%Y")) + try: + concDir = datahome + "CALPUFF_OUT/CALPUFF/" + year + "/m" + month + '/' + str(date.strftime("%Y%m%d")) + filenames, filePaths = concfiles(nConcFiles, concDir, SOX='SO2') + except AssertionError: + print(concDir + "does not exist skipping") + continue + for i, fname in enumerate(filePaths): + nhours += 1 # number of simulated hours + concA, concx = conc_array(ny, nx, fname, binLims) + concsum = concx + concsum + +# Average and write out +Concav = concsum/nhours +with open(year + month +'av.txt', 'w') as f: + for item in Concav: + f.write("%s\n" % item) diff --git a/Python/analysis_tools/monthlyP.py b/Python/analysis_tools/monthlyP.py new file mode 100644 index 0000000..270e42e --- /dev/null +++ b/Python/analysis_tools/monthlyP.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""pmaps +.. module:: genmaps + :platform: Unix + :synopis: +.. moduleauther: CEMAC (UoL) +.. description: This module was developed by CEMAC as part of the UNRESP + Project. This Script plots CALPUFF concrec data on a map. + :copyright: © 2019 University of Leeds. + :license: BSD-2 Clause. +Example: + To use:: + ./genemaps.py + - Date string, format YYYYMMDD, of the current CALPUFF run. + Used to locate directory containing the SO2 output files + (with assumed naming convention 'concrec0100**.dat', + where '**' goes from '01' through to '48') +.. CEMAC_UNRESPForcastingSystem: + https://github.com/cemac/UNRESPForcastingSystem +""" +from __future__ import print_function +import argparse +from dateutil.parser import parse +import datetime +from maptoolkit import * + + +def conc_array_plain(ny, nx, filePaths): + # Read in concentration data: + f = open(filePaths, 'r') + lines = f.read().splitlines() + f.close + # Process concentration data into desired format: + conc = np.array([float(X) for X in lines]) * 100**3 # ug/cm^3 -> ug/m^3 + concAry = np.reshape(conc, (ny, nx)) # Reshape data onto latlon grid + return concAry, conc + + +# READ IN COMMAND LINE ARGUMENTS +dstring = ("Date string, format YYYYMMDD, START date and End Date") +parser = argparse.ArgumentParser(description=dstring) +parser.add_argument("--start", help=dstring, type=str) +parser.add_argument("--end", help=dstring, type=str) +args = parser.parse_args() + +# Required files and info +start = datetime.datetime.strptime(args.start, "%Y%m%d") +end = datetime.datetime.strptime(args.end, "%Y%m%d") +datahome = "/nfs/see-fs-01_users/earhbu/UNRESP_SPACE/UNRESPForecastingSystem/" +nConcFiles = 48 +binLims = [10, 350, 600, 2600, 9000, 14000] # SO2 bin limits +xyFile = "../../data/xy_masaya.dat" +glat, glon, latMin, latMax, lonMin, lonMax, ny, nx = genxy(xyFile) +datelist = [start + datetime.timedelta(days=x) for x in range(0, (end-start).days)] +nhours = 0 +concsum = np.zeros(18000) +threshold = 100 +countersum = np.zeros((ny, nx)) +concAsum = np.zeros((ny, nx)) +concsum = np.zeros(18000) +concsumcount = np.zeros(18000) +for date in datelist: + month = str(date.strftime("%Y%m")) + year = str(date.strftime("%Y")) + try: + concDir = datahome + "CALPUFF_OUT/CALPUFF/" + year + "/m" + month + '/' + str(date.strftime("%Y%m%d")) + filenames, filePaths = concfiles(nConcFiles, concDir, SOX='SO2') + except AssertionError: + print(concDir + "does not exist skipping") + continue + for i, fname in enumerate(filePaths): + nhours += 1 # number of simulated hours + concA, concx = conc_array(ny, nx, fname, binLims) + concsum = concx + concsum + counter = np.zeros_like(concA) + counter[concA >= 100] = 1 + concAsum = concAsum + concA + countersum = counter + countersum + c = np.zeros(18000) + c[concx>=100] = 1 + concsumcount = c + concsumcount + +# Average and write out +Concav = concsumcount/nhours +with open(year + month +'avP.txt', 'w') as f: + for item in Concav: + f.write("%s\n" % item) diff --git a/Python/analysis_tools/monthly_stamps.py b/Python/analysis_tools/monthly_stamps.py new file mode 100644 index 0000000..9da648a --- /dev/null +++ b/Python/analysis_tools/monthly_stamps.py @@ -0,0 +1,106 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.image as mpimg +from mpl_toolkits.mplot3d import axes3d +import matplotlib as mpl +from mpl_toolkits.basemap import Basemap +from matplotlib.font_manager import FontProperties +import os +import datetime as dt +import pytz +import utm +import gmplot +from dateutil.parser import parse +from maptoolkit import * + + +class stamps(object): + def _init_(): + s.nConcFiles = 48 + s.binLims = [10, 350, 600, 2600, 9000, 14000] # SO2 bin limits + xyFile = "../../data/xy_masaya.dat" + s.glat, s.glon, s.latMin, s.latMax, s.lonMin, s.lonMax, s.ny, nx = genxy(xyFile) + im = gen_im(lonMin, latMin, lonMax, latMax, imtype="World_Shaded_Relief",) + topo = 'World_Shaded_Relief' + binLimsSO4 = [1E-8, 12, 35, 55, 150, 250] # SO4 bin limits from: + # http://mkwc.ifa.hawaii.edu/vmap/hysplit/ + colsHex = ['#FFFFFF', '#0cec0c', '#FFFF00', '#FF6600', '#FF0000', + '#800080', '#8F246B'] # Hex codes for SO2 colour bins + towns = (' El Panama', ' Rigoberto', ' Pacaya', ' El Crucero', + ' La Concepcion', ' Masaya', ' San Marcos', ' San Rafael del Sur', + ' Diriamba', ' Jinotepe', ' Masatepe') + townCoords = ((-86.2058, 11.972), (-86.2021, 11.9617), (-86.3013, 11.9553), + (-86.3113, 11.9923), (-86.189772, 11.936161), + (-86.096053, 11.973523), (-86.20317, 11.906584), + (-86.43639, 11.847034), (-86.239592, 11.85632), + (-86.19993, 11.85017), (-86.143758, 11.91512)) + cities = (' MANAGUA',) + cityCoords = ((-86.29, 12.12),) + volcCoords = (-86.1608, 11.9854) + font = FontProperties() + font.set_weight('bold') + font.set_family('monospace') + # SET BIN COLOURS + cmap = mpl.colors.ListedColormap(colsHex[1:-1]) + cmap.set_under(colsHex[0]) + cmap.set_over(colsHex[-1]) + normso4 = mpl.colors.BoundaryNorm(boundaries=binLimsSO4, ncolors=5) + norm = mpl.colors.BoundaryNorm(boundaries=binLims, ncolors=5) + glat, glon, latMin, latMax, lonMin, lonMax, ny, nx = genxy(xyFile) + fle = '2018201806av.txt' + + def plot_av_stamp(s, ita, im, fle, SOX): + """Plot static maps + """ + font = s.font + SOXf = r'SO$_' + SOX[-1] + '$' + so2title = ('Atmospheric ' + SOX + ' concentrations at ' + + 'ground level (hourly means). \n GCRF UNRESP') + plt.figure(figsize=(16, 12)) + if SOX == "SO4": + binLims = s.binLimsSO4 + norm = s.normso4 + else: + binLims = s.binLims + norm = s.norm + concA, concx = conc_array(s.ny, s.nx, fle, binLims) + latMin, latMax, lonMin = s.latMin, s.latMax, s.lonMin + lonMax = s.lonMax + bmap = Basemap(llcrnrlon=lonMin, llcrnrlat=latMin, + urcrnrlon=lonMax, urcrnrlat=latMax) + bmap.imshow(im, origin='upper') + bmap.pcolormesh(glon, glat, concA, + norm=norm, cmap=s.cmap, alpha=0.5) + cbar = bmap.colorbar(location='bottom', pad='20%', cmap=s.cmap, + norm=norm, boundaries=[0.] + binLims + + [100000.], extend='both', extendfrac='auto', + ticks=binLims, spacing='uniform') + cbar.ax.set_xticklabels(['v low', 'low', 'moderate', 'mod high', + 'high', 'v high']) # horizontal colorbar + cbar.set_label(label=(SOX + ' concentration'), fontsize=18) + cbar.ax.tick_params(labelsize=16) + cbar.solids.set(alpha=1) + latTicks = np.arange(round(latMin, 1), round(latMax, 1) + 0.1, 0.1) + lonTicks = np.arange(round(lonMin, 1), round(lonMax, 1) + 0.1, 0.2) + bmap.drawparallels(latTicks, labels=[1, 0, 0, 0], linewidth=0.0, + fontsize=16) + bmap.drawmeridians(lonTicks, labels=[0, 0, 0, 1], linewidth=0.0, + fontsize=16) + for i, town in enumerate(s.towns): + plt.plot(s.townCoords[i][0], s.townCoords[i] + [1], 'ok', markersize=4) + plt.text(s.townCoords[i][0], s.townCoords[i][1], town, + color='k', fontproperties=font, fontsize=12) + for i, city in enumerate(cities): + plt.plot(s.cityCoords[i][0], s.cityCoords[i] + [1], 'sk', markersize=6) + plt.text(s.cityCoords[i][0], s.cityCoords[i][1], city, + fontproperties=font, fontsize=16) + font0 = FontProperties() + font0.set_family('monospace') + plt.plot(s.volcCoords[0], volcCoords[1], '^r', markersize=6) + PNGfile = 'test.png' + print("Writing out file " + PNGfile) + PNGpath = os.path.join('images/', PNGfile) + plt.savefig(PNGpath, dpi=250) + plt.close() diff --git a/Python/analysis_tools/monthlygen.sh b/Python/analysis_tools/monthlygen.sh new file mode 100755 index 0000000..2ecf9a0 --- /dev/null +++ b/Python/analysis_tools/monthlygen.sh @@ -0,0 +1,31 @@ +#!/bin/bash - +#title :monthlygen.sh +#description :generated average conc files +#author :CEMAC - Helen +#date :20190429 +#version :1.0 +#usage :./monthlygen.sh +#notes : +#bash_version :4.2.46(2)-release +#============================================================================ + +python monthly.py --start 20180601 --end 20180630 +echo "June 2018" +python monthly.py --start 20180701 --end 20180731 +echo "July 2018" +python monthly.py --start 20180801 --end 20180831 +echo "August 2018" +python monthly.py --start 20180901 --end 20180930 +echo "Sept 2018" +python monthly.py --start 20181001 --end 20181031 +echo "Oct 2018" +python monthly.py --start 20181101 --end 20181130 +echo "Nov 2018" +python monthly.py --start 20181201 --end 20181231 +echo "Dec 2018" +python monthly.py --start 20190101 --end 20190131 +echo "Jan 2019" +python monthly.py --start 20190201 --end 20190228 +echo "Feb 2019" +python monthly.py --start 20190301 --end 20190331 +echo "March 2019" diff --git a/Python/analysis_tools/monthlygenP.sh b/Python/analysis_tools/monthlygenP.sh new file mode 100755 index 0000000..43e0079 --- /dev/null +++ b/Python/analysis_tools/monthlygenP.sh @@ -0,0 +1,31 @@ +#!/bin/bash - +#title :monthlygen.sh +#description :generated average conc files +#author :CEMAC - Helen +#date :20190429 +#version :1.0 +#usage :./monthlygen.sh +#notes : +#bash_version :4.2.46(2)-release +#============================================================================ + +python monthlyP.py --start 20180601 --end 20180630 +echo "June 2018" +python monthlyP.py --start 20180701 --end 20180731 +echo "July 2018" +python monthlyP.py --start 20180801 --end 20180831 +echo "August 2018" +python monthlyP.py --start 20180901 --end 20180930 +echo "Sept 2018" +python monthlyP.py --start 20181001 --end 20181031 +echo "Oct 2018" +python monthlyP.py --start 20181101 --end 20181130 +echo "Nov 2018" +python monthlyP.py --start 20181201 --end 20181231 +echo "Dec 2018" +python monthlyP.py --start 20190101 --end 20190131 +echo "Jan 2019" +python monthlyP.py --start 20190201 --end 20190228 +echo "Feb 2019" +python monthlyP.py --start 20190301 --end 20190331 +echo "March 2019" diff --git a/Python/analysis_tools/plotstamps.py b/Python/analysis_tools/plotstamps.py new file mode 100644 index 0000000..c30125c --- /dev/null +++ b/Python/analysis_tools/plotstamps.py @@ -0,0 +1,50 @@ +from __future__ import print_function +import argparse +import os +import matplotlib as mpl +from dateutil.parser import parse +import maptoolkit as mtk +# University System python may be broken +# If some one insists on using it... +backend = mpl.get_backend() +if backend == 'Qt4Agg' and sys.version_info[0] == 2: + # Fix the backend + print('swapping to Agg Backend') + print('Please consider using anaconda') + mpl.use('Agg') +# DO NOT MOVE ABOVE BACKEND FIX +import matplotlib.pyplot as plt # KEEP ME HERE!!! +################################ + +mpt = mtk.Masaya_Maps('20190412') +im = mtk.gen_im(mpt.lonMin, mpt.latMin, mpt.lonMax, + mpt.latMax, imtype='World_Shaded_Relief') +Figletter = ['a) ', 'b) ', 'c) ', 'd) ', 'e) ', 'f) ', 'g) ', 'h) ', 'j) ', + 'k) ', 'l) '] +filelist = ['2018201806av.txt', '2018201807av.txt', '2018201808av.txt', + '2018201809av.txt', '2018201810av.txt', '2018201811av.txt', + '2018201812av.txt', '2019201901av.txt', '2019201902av.txt', + '2019201903av.txt'] +mnth = ['Jun 2018', 'Jul 2018', 'Aug 2018', 'Sep 2018', 'Oct 2018', 'Nov 2018', + 'Dec 2018', 'Jan 2019', 'Feb 2019', 'Mar 2019'] +SOX = 'SO2' +fig = plt.figure(figsize=(30, 16)) +for i, fle in enumerate(filelist): + let = Figletter[i] + ax = fig.add_subplot(3, 4, i+1) + flepath = 'analysis_tools/' + fle + concA = mpt.plot_av_stamp(ax, mnth[i], let, im, flepath, SOX) +plt.tight_layout() +cax = fig.add_axes([1, 0.1, 0.03, 0.8]) +cbar = fig.colorbar(concA, cax=cax, cmap=mpt.cmap, norm=mpt.norm, boundaries=[0.] + mpt.binLims + + [100000.], extend='both', ticks=mpt.binLims) +cbar.ax.set_xticklabels(['v low', 'low', 'moderate', 'mod high', + 'high', 'v high']) # horizontal colorbar +cbar.set_label(label=(SOX + ' concentration'), fontsize=18) +cbar.ax.tick_params(labelsize=16) +cbar.solids.set(alpha=1) +PNGfile = 'Avconc_stamp.png' +print("Writing out file " + PNGfile) +PNGpath = os.path.join('analysis_tools/images/', PNGfile) +plt.savefig(PNGpath, dpi=250) +plt.close() diff --git a/Python/analysis_tools/pmap.py b/Python/analysis_tools/pmap.py new file mode 100755 index 0000000..98c9657 --- /dev/null +++ b/Python/analysis_tools/pmap.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""pmaps +.. module:: genmaps + :platform: Unix + :synopis: +.. moduleauther: CEMAC (UoL) +.. description: This module was developed by CEMAC as part of the UNRESP + Project. This Script plots CALPUFF concrec data on a map. + :copyright: © 2019 University of Leeds. + :license: BSD-2 Clause. +Example: + To use:: + ./genemaps.py + - Date string, format YYYYMMDD, of the current CALPUFF run. + Used to locate directory containing the SO2 output files + (with assumed naming convention 'concrec0100**.dat', + where '**' goes from '01' through to '48') +.. CEMAC_UNRESPForcastingSystem: + https://github.com/cemac/UNRESPForcastingSystem +""" +from __future__ import print_function +import argparse +from dateutil.parser import parse +import datetime +from maptoolkit import * +import netCDF4 + + +def conc_array_plain(ny, nx, filePaths): + # Read in concentration data: + f = open(filePaths, 'r') + lines = f.read().splitlines() + f.close + # Process concentration data into desired format: + conc = np.array([float(X) for X in lines]) * 100**3 # ug/cm^3 -> ug/m^3 + concAry = np.reshape(conc, (ny, nx)) # Reshape data onto latlon grid + return concAry, conc + + +# READ IN COMMAND LINE ARGUMENTS +dstring = ("Date string, format YYYYMMDD, START date and End Date") +parser = argparse.ArgumentParser(description=dstring) +parser.add_argument("--start", help=dstring, type=str) +parser.add_argument("--end", help=dstring, type=str) +args = parser.parse_args() + +# Required files and info +start = datetime.datetime.strptime(args.start, "%Y%m%d") +end = datetime.datetime.strptime(args.end, "%Y%m%d") +datahome = "/nfs/see-fs-01_users/earhbu/UNRESP_SPACE/UNRESPForecastingSystem/" +nConcFiles = 48 +binLims = [10, 350, 600, 2600, 9000, 14000] # SO2 bin limits +xyFile = "../data/xy_masaya.dat" +glat, glon, latMin, latMax, lonMin, lonMax, ny, nx = genxy(xyFile) +datelist = [start + datetime.timedelta(days=x) for x in range(0, (end-start).days)] +""" +the probability is simply calculated by checking each single hourly file and +"counting" how many time the concentration in each grid point is over a +specific threshold (here 100 micrograms/m3). This number (obtained for each +grid point) is then normalized over the total number of simulated hours.""" + +threshold = 100 + +nhours = 0 +countersum = np.zeros((ny, nx)) +concAsum = np.zeros((ny, nx)) +concsum = np.zeros(18000) +concsumcount = np.zeros(18000) +for date in datelist: + month = str(date.strftime("%Y%m")) + try: + concDir = datahome + "CALPUFF_OUT/CALPUFF/m" + month + '/' + str(date.strftime("%Y%m%d")) + filenames, filePaths = concfiles(nConcFiles, concDir, SOX='SO2') + except AssertionError: + print(concDir + "does not exist skipping") + continue + for i, fname in enumerate(filePaths): + nhours += 1 # number of simulated hours + concA, concx = conc_array(ny, nx, fname, binLims) + counter = np.zeros_like(concA) + counter[concA >= 100] = 1 + concAsum = concAsum + concA + countersum = counter + countersum + concsum = concx + concsum + c = np.zeros(18000) + c[concx>=100] = 1 + concsumcount = c + concsumcount +# Write Netcdf out +f = netCDF4.Dataset('conc.nc', 'w') +f.createDimension('lon', nx) +f.createDimension('lat', ny) +f.createDimension('N', 1) +f.createDimension('points', 18000) +h2 = f.createVariable('count', 'float', ('lon', 'lat')) +h2[:] = countersum +h1 = f.createVariable('conc', 'float', ('lon', 'lat')) +h1[:] = concAsum +h3 = f.createVariable('long', 'float', ('lon')) +h3[:] = glon[0, :] +h4 = f.createVariable('lati', 'float', ('lat')) +h4[:] = glat[:, 0] +h5 = f.createVariable('nhours', 'float', ('N')) +h5[:] = nhours +h6 = f.createVariable('concpoint', 'float', ('points')) +h6[:] = concsum +h7 = f.createVariable('countpoints', 'float', ('points')) +h7[:] = concsumcount +f.close() diff --git a/Python/analysis_tools/pmap_plot.py b/Python/analysis_tools/pmap_plot.py new file mode 100644 index 0000000..f491692 --- /dev/null +++ b/Python/analysis_tools/pmap_plot.py @@ -0,0 +1,85 @@ +import netCDF4 +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.image as mpimg +from mpl_toolkits.mplot3d import axes3d +import matplotlib as mpl +from mpl_toolkits.basemap import Basemap +from matplotlib.font_manager import FontProperties +import os +import datetime as dt +import pytz +import utm +import gmplot +import cartopy.crs as ccrs +import cartopy.io.img_tiles as cimgt +import cartopy +import cartopy.feature as cfeat +from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER +import netCDF4 +from dateutil.parser import parse +from maptoolkit import * +file2read = netCDF4.Dataset('conc.nc', 'r') +lat = file2read.variables['lati'][:] +lon = file2read.variables['long'][:] +nhours = file2read.variables['nhours'][:] +count = file2read.variables['countpoints'][:] +datahome = "/nfs/see-fs-01_users/earhbu/UNRESP_SPACE/UNRESPForecastingSystem/" +nConcFiles = 48 +binLims = [10, 350, 600, 2600, 9000, 14000] # SO2 bin limits +xyFile = "../data/xy_masaya.dat" +glat, glon, latMin, latMax, lonMin, lonMax, ny, nx = genxy(xyFile) +plt.figure(figsize=(16, 12)) +im = gen_im(lonMin, latMin, lonMax, latMax, imtype="World_Shaded_Relief",) +so2title = (r'june 2018 - feb 2019. probability greater than 100ug/m$^3$'+' \n (cut off reduced to 0.001)') +bmap = Basemap(llcrnrlon=lonMin, llcrnrlat=latMin, + urcrnrlon=lonMax, urcrnrlat=latMax) +bmap.imshow(im, origin='upper') +concAry = np.reshape(count, (ny, nx)) +#concAry[concAry < 0.05] = np.nan +concmap = concAry / nhours*100 +concmap = np.ma.masked_array(concmap, concmap < 0.001) +bmap.contourf(glon, glat, concmap, cmap=plt.cm.plasma_r) +cbar = bmap.colorbar(location='bottom', pad='20%', extend='both', extendfrac='auto', cmap=plt.cm.plasma_r) +cbar.set_label(label=('probability %'), fontsize=18) +cbar.ax.tick_params(labelsize=16) +cbar.solids.set(alpha=1) +latTicks = np.arange(round(latMin, 1), round(latMax, 1) + 0.1, 0.1) +lonTicks = np.arange(round(lonMin, 1), round(lonMax, 1) + 0.1, 0.2) +bmap.drawparallels(latTicks, labels=[1, 0, 0, 0], linewidth=0.0, + fontsize=16) +bmap.drawmeridians(lonTicks, labels=[0, 0, 0, 1], linewidth=0.0, + fontsize=16) +towns = (' El Panama', ' Rigoberto', ' Pacaya', ' El Crucero', + ' La Concepcion', ' Masaya', ' San Marcos', + ' San Rafael del Sur', ' Diriamba', ' Jinotepe', + ' Masatepe') +townCoords = ((-86.2058, 11.972), (-86.2021, 11.9617), + (-86.3013, 11.9553), (-86.3113, 11.9923), + (-86.189772, 11.936161), (-86.096053, 11.973523), + (-86.20317, 11.906584), (-86.43639, 11.847034), + (-86.239592, 11.85632), (-86.19993, 11.85017), + (-86.143758, 11.91512)) +cities = (' MANAGUA',) +cityCoords = ((-86.29, 12.12),) +volcCoords = (-86.1608, 11.9854) +font = FontProperties() +font.set_weight('bold') +font.set_family('monospace') +tc = 'k' +for i, town in enumerate(towns): + plt.plot(townCoords[i][0], townCoords[i] + [1], 'ok', markersize=4) + plt.text(townCoords[i][0], townCoords[i][1], town, + color=tc, fontproperties=font, fontsize=12) +for i, city in enumerate(cities): + plt.plot(cityCoords[i][0], cityCoords[i] + [1], 'sk', markersize=6) + plt.text(cityCoords[i][0], cityCoords[i][1], city, + fontproperties=font, fontsize=16) +font0 = FontProperties() +font0.set_family('monospace') +plt.plot(volcCoords[0], volcCoords[1], '^r', markersize=6) +plt.suptitle(so2title, fontsize=24) +plt.savefig('test.png', dpi=250) +plt.close() diff --git a/Python/maptoolkit.py b/Python/maptoolkit.py index e5f0336..a9798f9 100644 --- a/Python/maptoolkit.py +++ b/Python/maptoolkit.py @@ -358,6 +358,92 @@ def plot_googlemap1(s, ita, filePaths, SOX): print("Writing out file " + HTMLfile) gmap.draw(os.path.join(s.outDir, HTMLfile)) + def plot_av_stamp(s, ax, mnth, let, im, fle, SOX): + """Plot static maps + """ + font = s.font + SOXf = r'SO$_' + SOX[-1] + '$' + so2title = ('Atmospheric ' + SOX + ' concentrations at ' + + 'ground level (hourly means). \n GCRF UNRESP') + if SOX == "SO4": + binLims = s.binLimsSO4 + norm = s.normso4 + else: + binLims = s.binLims + norm = s.norm + concA, concx = conc_array(s.ny, s.nx, fle, binLims) + latMin, latMax, lonMin = s.latMin, s.latMax, s.lonMin + lonMax = s.lonMax + ax.bmap = Basemap(llcrnrlon=lonMin, llcrnrlat=latMin, + urcrnrlon=lonMax, urcrnrlat=latMax) + ax.bmap.imshow(im, origin='upper') + p = ax.bmap.pcolormesh(s.glon, s.glat, concA / 100**3, + norm=norm, cmap=s.cmap, alpha=0.5) + latTicks = np.arange(round(latMin, 1), round(latMax, 1) + 0.1, 0.1) + lonTicks = np.arange(round(lonMin, 1), round(lonMax, 1) + 0.1, 0.2) + ax.bmap.drawparallels(latTicks, labels=[1, 0, 0, 0], linewidth=0.0, + fontsize=16) + ax.bmap.drawmeridians(lonTicks, labels=[0, 0, 0, 1], linewidth=0.0, + fontsize=16) + for i, town in enumerate(s.towns): + ax.plot(s.townCoords[i][0], s.townCoords[i] + [1], 'ok', markersize=4) + ax.text(s.townCoords[i][0], s.townCoords[i][1], town, + color='k', fontproperties=font, fontsize=8) + for i, city in enumerate(s.cities): + ax.plot(s.cityCoords[i][0], s.cityCoords[i] + [1], 'sk', markersize=6) + ax.text(s.cityCoords[i][0], s.cityCoords[i][1], city, + fontproperties=font, fontsize=16) + font0 = FontProperties() + font0.set_family('monospace') + ax.plot(s.volcCoords[0], s.volcCoords[1], '^r', markersize=6) + ax.set_title(str(let)+str(mnth), fontsize=20) + return p + + def plot_av_Prob(s, ax, mnth, let, im, fle, SOX): + """Plot static maps + """ + font = s.font + SOXf = r'SO$_' + SOX[-1] + '$' + so2title = ('Atmospheric ' + SOX + ' concentrations at ' + + 'ground level (hourly means). \n GCRF UNRESP') + if SOX == "SO4": + binLims = s.binLimsSO4 + norm = s.normso4 + else: + binLims = s.binLims + norm = s.norm + binLims[0] = 0 + concA, concx = conc_array(s.ny, s.nx, fle, binLims) + latMin, latMax, lonMin = s.latMin, s.latMax, s.lonMin + lonMax = s.lonMax + ax.bmap = Basemap(llcrnrlon=lonMin, llcrnrlat=latMin, + urcrnrlon=lonMax, urcrnrlat=latMax) + ax.bmap.imshow(im, origin='upper') + p = ax.bmap.pcolormesh(s.glon, s.glat, concA * 1000, cmap=plt.cm.plasma_r) + latTicks = np.arange(round(latMin, 1), round(latMax, 1) + 0.1, 0.1) + lonTicks = np.arange(round(lonMin, 1), round(lonMax, 1) + 0.1, 0.2) + ax.bmap.drawparallels(latTicks, labels=[1, 0, 0, 0], linewidth=0.0, + fontsize=16) + ax.bmap.drawmeridians(lonTicks, labels=[0, 0, 0, 1], linewidth=0.0, + fontsize=16) + for i, town in enumerate(s.towns): + ax.plot(s.townCoords[i][0], s.townCoords[i] + [1], 'ok', markersize=4) + ax.text(s.townCoords[i][0], s.townCoords[i][1], town, + color='k', fontproperties=font, fontsize=8) + for i, city in enumerate(s.cities): + ax.plot(s.cityCoords[i][0], s.cityCoords[i] + [1], 'sk', markersize=6) + ax.text(s.cityCoords[i][0], s.cityCoords[i][1], city, + fontproperties=font, fontsize=16) + font0 = FontProperties() + font0.set_family('monospace') + ax.plot(s.volcCoords[0], s.volcCoords[1], '^r', markersize=6) + ax.set_title(str(let)+str(mnth), fontsize=20) + return p + def plot_diff(s): """ """ diff --git a/Python/plotstamps.py b/Python/plotstamps.py new file mode 100644 index 0000000..c30125c --- /dev/null +++ b/Python/plotstamps.py @@ -0,0 +1,50 @@ +from __future__ import print_function +import argparse +import os +import matplotlib as mpl +from dateutil.parser import parse +import maptoolkit as mtk +# University System python may be broken +# If some one insists on using it... +backend = mpl.get_backend() +if backend == 'Qt4Agg' and sys.version_info[0] == 2: + # Fix the backend + print('swapping to Agg Backend') + print('Please consider using anaconda') + mpl.use('Agg') +# DO NOT MOVE ABOVE BACKEND FIX +import matplotlib.pyplot as plt # KEEP ME HERE!!! +################################ + +mpt = mtk.Masaya_Maps('20190412') +im = mtk.gen_im(mpt.lonMin, mpt.latMin, mpt.lonMax, + mpt.latMax, imtype='World_Shaded_Relief') +Figletter = ['a) ', 'b) ', 'c) ', 'd) ', 'e) ', 'f) ', 'g) ', 'h) ', 'j) ', + 'k) ', 'l) '] +filelist = ['2018201806av.txt', '2018201807av.txt', '2018201808av.txt', + '2018201809av.txt', '2018201810av.txt', '2018201811av.txt', + '2018201812av.txt', '2019201901av.txt', '2019201902av.txt', + '2019201903av.txt'] +mnth = ['Jun 2018', 'Jul 2018', 'Aug 2018', 'Sep 2018', 'Oct 2018', 'Nov 2018', + 'Dec 2018', 'Jan 2019', 'Feb 2019', 'Mar 2019'] +SOX = 'SO2' +fig = plt.figure(figsize=(30, 16)) +for i, fle in enumerate(filelist): + let = Figletter[i] + ax = fig.add_subplot(3, 4, i+1) + flepath = 'analysis_tools/' + fle + concA = mpt.plot_av_stamp(ax, mnth[i], let, im, flepath, SOX) +plt.tight_layout() +cax = fig.add_axes([1, 0.1, 0.03, 0.8]) +cbar = fig.colorbar(concA, cax=cax, cmap=mpt.cmap, norm=mpt.norm, boundaries=[0.] + mpt.binLims + + [100000.], extend='both', ticks=mpt.binLims) +cbar.ax.set_xticklabels(['v low', 'low', 'moderate', 'mod high', + 'high', 'v high']) # horizontal colorbar +cbar.set_label(label=(SOX + ' concentration'), fontsize=18) +cbar.ax.tick_params(labelsize=16) +cbar.solids.set(alpha=1) +PNGfile = 'Avconc_stamp.png' +print("Writing out file " + PNGfile) +PNGpath = os.path.join('analysis_tools/images/', PNGfile) +plt.savefig(PNGpath, dpi=250) +plt.close() diff --git a/Python/plotstampsP.py b/Python/plotstampsP.py new file mode 100644 index 0000000..908cb08 --- /dev/null +++ b/Python/plotstampsP.py @@ -0,0 +1,42 @@ +from __future__ import print_function +import argparse +import os +import matplotlib as mpl +from dateutil.parser import parse +import maptoolkit as mtk +# University System python may be broken +# If some one insists on using it... +backend = mpl.get_backend() +if backend == 'Qt4Agg' and sys.version_info[0] == 2: + # Fix the backend + print('swapping to Agg Backend') + print('Please consider using anaconda') + mpl.use('Agg') +# DO NOT MOVE ABOVE BACKEND FIX +import matplotlib.pyplot as plt # KEEP ME HERE!!! +################################ + +mpt = mtk.Masaya_Maps('20190412') +im = mtk.gen_im(mpt.lonMin, mpt.latMin, mpt.lonMax, + mpt.latMax, imtype='World_Shaded_Relief') +Figletter = ['a) ', 'b) ', 'c) ', 'd) ', 'e) ', 'f) ', 'g) ', 'h) ', 'j) ', + 'k) ', 'l) '] +filelist = ['2018201806avP.txt', '2018201807avP.txt', '2018201808avP.txt', + '2018201809avP.txt', '2018201810avP.txt', '2018201811avP.txt', + '2018201812avP.txt', '2019201901avP.txt', '2019201902avP.txt', + '2019201903avP.txt'] +mnth = ['Jun 2018', 'Jul 2018', 'Aug 2018', 'Sep 2018', 'Oct 2018', 'Nov 2018', + 'Dec 2018', 'Jan 2019', 'Feb 2019', 'Mar 2019'] +SOX = 'SO2' +fig = plt.figure(figsize=(30, 16)) +for i, fle in enumerate(filelist): + let = Figletter[i] + ax = fig.add_subplot(3, 4, i+1) + flepath = 'analysis_tools/' + fle + concA = mpt.plot_av_stamp(ax, mnth[i], let, im, flepath, SOX) +plt.tight_layout() +PNGfile = 'P_stamp.png' +print("Writing out file " + PNGfile) +PNGpath = os.path.join('analysis_tools/images/', PNGfile) +plt.savefig(PNGpath, dpi=250) +plt.close()