io.py 4.25 KB
Newer Older
Eva Börgens's avatar
Eva Börgens committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
import netCDF4
import re
import numpy as np
import scipy.special as sp
import math
import sys

def read_netcdf(filename):
    try:
        with netCDF4.Dataset(filename, 'r') as S:
            lon = S.variables['lon'][:]
            lat = S.variables['lat'][:]
            time = S.variables['time'][:]
            grid = S.variables['tws'][:]
            grid_std = S.variables['std_tws'][:]
        grid_std = np.array([np.nanmean(g) for g in grid_std])
    except Exception:
        print(" Failed to read %s. Variable names in file have to be: lon, lat, time, tws, std_tws" %filename)
        sys.exit()

    return grid, grid_std, lon, lat, time

def read_ascii(filename):
    try:
        with open(filename) as fid:
            coords = []
            for l in fid:
                if l[0]=='#':
                    continue
                line = [float(i) for i in re.split('\s+', l.strip())]
                coords.append([line[0], line[1]])
    except:
        print("Could not read %s: Please ensure that every line contains two columns and comment lines are indicated with '#'" %filename)
        sys.exit()
    return coords

def test_coordiantes(filename):
    coords = read_ascii(filename)

    lon = np.array([r[0] for r in coords])
    lat = np.array([r[1] for r in coords])
    if len(lon)<2:
        raise RuntimeError("Coordinate list of region have to contain at least 2 points")

    u_lon = np.unique(lon)
    u_lat = np.unique(lat)

    if not same_dist_elems(u_lon):
        raise RuntimeError("Longitudes of region not evenly spaced")
    if not same_dist_elems(u_lat):
        raise RuntimeError("Latitudes of region not evenly spaced")

    return

def same_dist_elems(arr):
    diff = arr[1] - arr[0]
    for x in range(1, len(arr) - 1):
        if arr[x + 1] - arr[x] != diff:
            return False
    return True

def save_results(outname, inname, results, region_coords, timeseries, flag_uncertainty, flag_matrix, flag_timeseries):
    lon = np.array([r[0] for r in region_coords])
    lat = np.array([r[1] for r in region_coords])

    with netCDF4.Dataset(inname, format='NETCDF4') as Sin:
        time = Sin.variables['time'][:]
        time_unit = Sin.variables['time'].units
        time_calendar = Sin.variables['time'].calendar

        tws_unit = Sin.variables['tws'].units

    try:
        with netCDF4.Dataset(outname, 'w', format='NETCDF4') as Sout:
            Sout.createDimension('time', len(time))
            time_out = Sout.createVariable('time', 'i4', 'time')
            time_out.units = time_unit
            time_out.calendar = time_calendar
            time_out[:] = time

            Sout.createDimension('pointID', len(region_coords))
            coords_id = Sout.createVariable('pointID', 'i4', 'pointID')
            coords_id[:] = np.arange(0, len(region_coords))

            coords_lon = Sout.createVariable('lon', 'f4', ('pointID'))
            coords_lon[:] = lon
            coords_lon.long_name = "Lon coordinates tuple of the grid points of the region"
            coords_lon.units = "degrees_north"

            coords_lat = Sout.createVariable('lat', 'f4', ('pointID'))
            coords_lat[:] = lat
            coords_lat.long_name = "Lat coordinates tuple of the grid points of the region"
            coords_lat.units = "degrees_east"

            if flag_uncertainty:
                uncertainty = Sout.createVariable('std', 'f4', 'time')
                uncertainty.units = tws_unit
                uncertainty.long_name = "standard deviation of regions mean tws time series"
                uncertainty[:] = results['uncertainty']
            if flag_matrix:
                cov_matrix = Sout.createVariable('cov_matrix', 'f4', ('time','pointID', 'pointID'))
                cov_matrix.long_name = "Covariance matrix between all points of the region"
                cov_matrix.units = tws_unit + '*' + tws_unit
                cov_matrix[:] = results['matrix']
            if flag_timeseries:
                timeseries_out = Sout.createVariable('tws_timeseries', 'f4', 'time')
                timeseries_out.units = tws_unit
                timeseries_out.long_name = "Regions mean tws time series"
                timeseries_out[:] = timeseries
    except PermissionError:
        print("Could not open %s for writing results" %outname)
        sys.exit()
    return