diff --git a/.gitignore b/.gitignore index 669169d..99018d0 100644 --- a/.gitignore +++ b/.gitignore @@ -11,8 +11,9 @@ doc/build/* .project .pydev* f2py_build/f90/* -FortFlex.so +FortFlex*.so .idea/* /.cache/ .*.swp + diff --git a/reflexible/conv2netcdf4/flexpart_read.py b/reflexible/conv2netcdf4/flexpart_read.py index 6b6ea4c..7c0eadd 100644 --- a/reflexible/conv2netcdf4/flexpart_read.py +++ b/reflexible/conv2netcdf4/flexpart_read.py @@ -5,6 +5,7 @@ import os import datetime import re +import string import numpy as np import bcolz from math import pi, sqrt, cos @@ -40,67 +41,121 @@ def read_releases_v10(pathname): Returns ------- - A ctable object from bcolz package. + A dictionary with the release data + + .. tabularcolumns:: |l|L| + + =================== ======================================== + Keys Description + =================== ======================================== + releases release data object from bcolz package + spec_ids array of species id integers + =================== ======================================== """ - # Setup the container for the data - dtype = [('IDATE1', np.int32), ('ITIME1', np.int32), - ('IDATE2', np.int32), ('ITIME2', np.int32), - ('LON1', np.float32), ('LON2', np.float32), - ('LAT1', np.float32), ('LAT2', np.float32), - ('Z1', np.float32), ('Z2', np.float32), - ('ZKIND', np.int8), ('MASS', np.float32), - ('PARTS', np.int32), ('COMMENT', 'S32')] + # Read the input + input_str = open(pathname, 'r').read() + group_names = ['&RELEASES_CTRL\n', '&RELEASE\n'] + end_group = '/\n' + dtypes = [[('NSPEC', np.int32), ('SPECNUM_REL', object)], + [('IDATE1', np.int32), ('ITIME1', np.int32), + ('IDATE2', np.int32), ('ITIME2', np.int32), + ('LON1', np.float32), ('LON2', np.float32), + ('LAT1', np.float32), ('LAT2', np.float32), + ('Z1', np.float32), ('Z2', np.float32), + ('ZKIND', np.int8), ('MASS', np.float32), + ('PARTS', np.int32), ('COMMENT', 'S32')]] + + # Matches a key = value pair to obtain the value + namelist_group_re = r'[^\s,\=]+?\s*=\s*((?:[^=]+)|(?:(?:\"|\').+))(?=[\s,][^\s,\=]+\s*\=|$)' + cparams = bcolz.cparams(cname="lz4", clevel=6, shuffle=1) - ctable = bcolz.zeros(0, dtype=dtype, cparams=cparams) - nrecords = ctable['IDATE1'].chunklen - releases = np.zeros(nrecords, dtype=dtype) + + for k, group_name in enumerate(group_names): + len_marker = len(group_name) + ctable = bcolz.zeros(0, dtype=dtypes[k], cparams=cparams) + nrecords = ctable[ctable.names[0]].chunklen + values = np.zeros(nrecords, dtype=dtypes[k]) + if k == 0: + releases_species = ctable + elif k == 1: + releases = ctable + # Loop over all the marker groups + i, n = 0, 0 + while True: + i = input_str.find(group_name, i) + if i == -1: + break # group marker is not found anymore + j = input_str.find(end_group, i + 1) + n += 1 + group_block = input_str[i + len_marker: j] + i = j + values_tuple = tuple([s.strip(',"\''+string.whitespace) for s in re.findall(namelist_group_re, group_block)]) - # Prepare for reading the input - input_str = open(pathname, 'r').read() - marker = "&RELEASE\n" - len_marker = len(marker) - release_re = r'\S+=\s+[\"|\s](\S+)[,|\"|\w]' - - # Loop over all the marker groups - i, n = 0, 0 - while True: - i = input_str.find(marker, i) - j = input_str.find(marker, i + 1) - n += 1 - group_block = input_str[i + len_marker: j] - i = j - values = tuple(re.findall(release_re, group_block)) - try: - releases[(n - 1) % nrecords] = values - except ValueError: - print("Problem at: group: %d, %s" % (n, group_block)) - print("values:", values) - raise - if (n % nrecords) == 0: - ctable.append(releases) - if (i == -1) or (j == -1): - break # marker is not found anymore - # Remainder - ctable.append(releases[:n % nrecords]) - ctable.flush() + try: + values[(n - 1) % nrecords] = values_tuple + except ValueError: + print("Problem at: group: %d, %s" % (n, group_block)) + print("values:", values_tuple) + raise + if n % nrecords == 0: + ctable.append(values) - return ctable + # Remainder + ctable.append(values[:n % nrecords]) + ctable.flush() + + return {'releases': releases, + 'spec_ids': np.array(releases_species['SPECNUM_REL'][0].split(','), dtype=int)} -def read_releases_v9(path): +def read_releases_v9(path, headerrows=11): """Read metadata from a RELEASES path and return it as a dict. - Only 'release_point_names' entry returned. + Parameters + ---------- + path : pathname + Release file name + + Returns + ------- + A dictionary with the release metadata + + .. tabularcolumns:: |l|L| + + =================== ======================================== + Keys Description + =================== ======================================== + release_point_names array of release point name strings + spec_ids array of species id integers + =================== ======================================== """ + spec_start_line = headerrows+1 + old_format = False rpnames = [] with open(path) as f: prev_line = None + lnum = 0 + nspec = 0 + spec_ids = list() for line in f: - if prev_line is not None and "comment" in line: - rpnames.append(prev_line.strip()) - prev_line = line - # Return just the release point names for now - return {"release_point_names": np.array(rpnames, dtype="S45")} + lnum += 1 + if lnum == spec_start_line: + nspec = int(line.split()[0]) + elif lnum == spec_start_line+1 and 'Total' in line: + old_format = True + skip_lines = 2 + if not old_format: + if lnum > spec_start_line and lnum <= spec_start_line + nspec: + spec_ids.append(int(line.split()[0])) + else: + if prev_line is not None and "comment" in line: + rpnames.append(prev_line.strip()) + shift_lnum = lnum - spec_start_line - skip_lines + if shift_lnum > 0 and shift_lnum <= nspec*(skip_lines+1) and shift_lnum % (skip_lines+1) == 1: + spec_ids.append(int(line.split()[0])) + prev_line = line + return {'release_point_names': np.array(rpnames, dtype='S45'), + 'spec_ids': np.array(spec_ids, dtype=int)} def read_releases_old(path, headerrows=11): diff --git a/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh b/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh index 9f3ce4e..3c4d821 100755 --- a/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh +++ b/reflexible/conv2netcdf4/fortflex/build_FortFlex.sh @@ -1,7 +1,7 @@ #!/bin/bash -rm -f FortFlex.pyf f2py -m FortFlex -h FortFlex.pyf FortFlex.f f2py -c --fcompiler=gfortran FortFlex.pyf FortFlex.f +rm -f FortFlex.pyf mv *.so .. diff --git a/reflexible/conv2netcdf4/fortflex/build_FortFlex_py3.sh b/reflexible/conv2netcdf4/fortflex/build_FortFlex_py3.sh new file mode 100755 index 0000000..c6a95b3 --- /dev/null +++ b/reflexible/conv2netcdf4/fortflex/build_FortFlex_py3.sh @@ -0,0 +1,7 @@ +#!/bin/bash + +f2py3 -m FortFlex -h FortFlex.pyf FortFlex.f +f2py3 -c --fcompiler=gfortran FortFlex.pyf FortFlex.f +rm -f FortFlex.pyf +mv *.so .. + diff --git a/reflexible/conv2netcdf4/grid_read.py b/reflexible/conv2netcdf4/grid_read.py index 5a0e8c8..0b9496e 100644 --- a/reflexible/conv2netcdf4/grid_read.py +++ b/reflexible/conv2netcdf4/grid_read.py @@ -613,7 +613,7 @@ def readgridV8(H, **kwargs): # What species to return? nspec_ret = OPS.nspec_ret - if isinstance(nspec_ret, int): + if isinstance(nspec_ret, (int, np.int64)): nspec_ret = [nspec_ret] assert iter(nspec_ret), "nspec_ret must be iterable." @@ -622,7 +622,7 @@ def readgridV8(H, **kwargs): if OPS.time_ret is not None: get_dates = [] time_ret = OPS.time_ret - if isinstance(time_ret, int) == True: + if isinstance(time_ret, (int, np.int64)): time_ret = [time_ret] if time_ret[0] < 0: @@ -861,7 +861,7 @@ def readgridV6(H, **kwargs): # What species to return? nspec_ret = OPS.nspec_ret - if isinstance(nspec_ret, int): + if isinstance(nspec_ret, (int, np.int64)): nspec_ret = [nspec_ret] assert iter(nspec_ret), "nspec_ret must be iterable." @@ -870,7 +870,7 @@ def readgridV6(H, **kwargs): if OPS.time_ret is not None: get_dates = [] time_ret = OPS.time_ret - if isinstance(time_ret, int) == True: + if isinstance(time_ret, (int, np.int64)): time_ret = [time_ret] if time_ret[0] < 0: @@ -1097,7 +1097,7 @@ def fill_grids(H, nspec=0, FD=None): # assert H.direction == 'backward', "fill_backward is only valid for backward runs" # make sure npsec is iterable - if isinstance(nspec, int): + if isinstance(nspec, (int, np.int64)): species = [nspec] else: species = nspec diff --git a/reflexible/conv2netcdf4/tests/test_netcdf4_structure.py b/reflexible/conv2netcdf4/tests/test_netcdf4_structure.py index bf488d7..c769291 100644 --- a/reflexible/conv2netcdf4/tests/test_netcdf4_structure.py +++ b/reflexible/conv2netcdf4/tests/test_netcdf4_structure.py @@ -179,7 +179,7 @@ def test_LAGE(self): def test_species_mr(self): attr_names = ('units', 'long_name', 'decay', 'weightmolar', 'ohreact', 'kao', 'vsetaver', 'spec_ass') - for i in range(0, self.H.nspec): + for i in range(self.H.nspec): anspec = "%3.3d" % (i + 1) if self.ncid.iout in (1, 3, 5): var_name = "spec" + anspec + "_mr" @@ -194,7 +194,7 @@ def test_species_mr(self): def test_species_pptv(self): attr_names = ('units', 'long_name', 'decay', 'weightmolar', 'ohreact', 'kao', 'vsetaver', 'spec_ass') - for i in range(0, self.H.nspec): + for i in range(self.H.nspec): anspec = "%3.3d" % (i + 1) if self.ncid.iout in (2, 3): var_name = "spec" + anspec + "_pptv" @@ -217,7 +217,7 @@ def setup(self, request, tmpdir): def test_WDspecies(self): attr_names = ('units', 'weta', 'wetb', 'weta_in', 'wetb_in', 'wetc_in', 'wetd_in', 'dquer', 'henry') - for i in range(0, self.H.nspec): + for i in range(self.H.nspec): anspec = "%3.3d" % (i + 1) if self.dataset.wetdep: var_name = "WD_spec" + anspec @@ -232,7 +232,7 @@ def test_WDspecies(self): def test_DDspecies(self): attr_names = ('units', 'dryvel', 'reldiff', 'henry', 'f0', 'dquer', 'density', 'dsigma') - for i in range(0, self.H.nspec): + for i in range(self.H.nspec): anspec = "%3.3d" % (i + 1) if self.dataset.drydep: var_name = "DD_spec" + anspec diff --git a/reflexible/flexpart.py b/reflexible/flexpart.py index b6801d0..9fa93c5 100644 --- a/reflexible/flexpart.py +++ b/reflexible/flexpart.py @@ -48,7 +48,7 @@ def __init__(self, pathnames, nested=False): # Other config files self.Command = read_conffiles("COMMAND", self.fp_options, None) self.Releases = read_conffiles("RELEASES", self.fp_options, None) - self.Species = read_species(self.fp_options, self.Header.nspec) + self.Species = read_species(self.fp_options, self.Header.nspec, spec_ids=self.Releases['spec_ids']) def __str__(self): return "Flexpart('%s', nested=%s)" % (self.pathnames, self.nested) diff --git a/reflexible/scripts/create_ncfile.py b/reflexible/scripts/create_ncfile.py index 48c15ed..d6b4709 100755 --- a/reflexible/scripts/create_ncfile.py +++ b/reflexible/scripts/create_ncfile.py @@ -41,7 +41,7 @@ MIN_SIZE = False -def read_species(options_dir, nspec): +def read_species(options_dir, nspec, spec_ids=None): """Read metadata from SPECIES dir and return it as a dict. TODO: maybe the current version is specific of FLEXPART 9. @@ -72,7 +72,9 @@ def read_species(options_dir, nspec): } dspec = defaultdict(list) - for spec in range(1, nspec+1): + if spec_ids is None: + spec_ids = range(1, nspec+1) + for spec in spec_ids: path = os.path.join(species_dir, "SPECIES_%03d" % spec) with open(path) as fspec: for line in fspec: @@ -387,7 +389,7 @@ def write_header(H, ncid, wetdep, drydep, write_releases, species): chunksizes = (H.numxgrid, H.numygrid, H.numzgrid, 1, 1, 1)[::-1] dep_chunksizes = (H.numxgrid, H.numygrid, 1, 1, 1)[::-1] - for i in range(0, H.nspec): + for i in range(H.nspec): anspec = "%3.3d" % (i + 1) # iout: 1 conc. output (ng/m3), 2 mixing ratio (pptv), 3 both, # 4 plume traject, 5=1+4 @@ -608,7 +610,7 @@ def create_ncfile(pathnames, nested, wetdep=False, drydep=False, if options_dir: command = read_conffiles("COMMAND", options_dir, command_path) releases = read_conffiles("RELEASES", options_dir, releases_path) - species = read_species(options_dir, H.nspec) + species = read_species(options_dir, H.nspec, spec_ids=releases['spec_ids']) if outfile: # outfile has priority over previous flags diff --git a/reflexible/tests/test_flexpart.py b/reflexible/tests/test_flexpart.py index 4be411e..2189075 100644 --- a/reflexible/tests/test_flexpart.py +++ b/reflexible/tests/test_flexpart.py @@ -35,15 +35,15 @@ def test_Header(self): def test_releases(self): # Only test a few FP runs here if self.dataset.fp_name == 'Fwd1_V10.1': - releases = self.fprun.Releases[:] + releases = self.fprun.Releases['releases'] assert len(releases) == 1 assert tuple(releases[0]) == ( 20050501, 0, 20050502, 0, 10.5, 11.0, 59.75, 60.25, 0.0, 100.0, 1, 1.0, 10000, b'OsloRelease') elif self.dataset.fp_name == 'Fwd1_V9.02': - releases = self.fprun.Releases + releases = self.fprun.Releases['release_point_names'] assert len(releases) == 1 - assert releases['release_point_names'][0] == b'RELEASE_TEST1' + assert releases[0] == b'RELEASE_TEST1' def test_command(self): command = self.fprun.Command @@ -73,4 +73,56 @@ def test_command(self): def test_species(self): species = self.fprun.Species - assert species['dryvel'] == [-9.99] + if self.dataset.fp_name == 'Fwd1_V10.1': + # species id = 15 + assert species == { + 'decay': [691200.0], + # TODO: Search string in reflexible.scripts.create_ncfile.read_species + # needs to be updated to handle below/in-cloud scavenging in V10 + #'weta': [1.0E-04], + #'wetb': [0.80], + 'reldiff': [-9.9], + 'henry': [-9.9], + 'f0': [-9.9], + 'dquer': [6.0E-7], + 'dsigma': [3.0E-1], + 'dryvel': [-9.99], + 'weightmolar': [-9.99], + 'ohreact': [-9.9E-09], + 'spec_ass': [-9], + 'kao': [-99.99] + } + elif self.dataset.fp_name == 'Fwd1_V9.02': + # species id = 1 + assert species == { + 'decay': [-999.9], + 'weta': [-9.9E-09], + 'wetb': [-9.9E-09], + 'reldiff': [-9.9], + 'henry': [-9.9], + 'f0': [-9.9], + 'dquer': [-9.9], + 'dsigma': [-9.9], + 'dryvel': [-9.99], + 'weightmolar': [350.00], + 'ohreact': [-9.9E-09], + 'spec_ass': [-9], + 'kao': [-99.99] + } + elif self.dataset.fp_name == 'Fwd2_V9.02': + # species id = 12 + assert species == { + 'decay': [-999.9], + 'weta': [5.0E-06], + 'wetb': [0.62], + 'reldiff': [-9.9], + 'henry': [1.0E-09], + 'f0': [1.0E-09], + 'dquer': [4.0E-7], + 'dsigma': [3.0E-1], + 'dryvel': [-9.99], + 'weightmolar': [-9.99], + 'ohreact': [-9.9E-09], + 'spec_ass': [-9], + 'kao': [-99.99] + }