__init__.py 23.3 KB
Newer Older
CEVAER's avatar
CEVAER committed
1
#!/usr/bin/env python
2
# -*- coding: utf-8 -*-
3
4
5
6
from __future__ import division
from builtins import str
from builtins import range
from past.utils import old_div
CEVAER's avatar
CEVAER committed
7
8
9
10
11
12
13
14
import sys
import os
import datetime
import glob
import netCDF4
import numpy as np
import numpy.ma as ma
import logging
ARCHER's avatar
ARCHER committed
15
import warnings
16
from collections import OrderedDict
CEVAER's avatar
CEVAER committed
17
18
19
20
21
22

logging.basicConfig()
logger = logging.getLogger(__name__)
logger.setLevel(logging.DEBUG)

def readFile(fname):
ARCHER's avatar
ARCHER committed
23
24
25
26
27
28
29
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        data={}
        #list_var=('owiLon','owiLat','owiHeading','owiEcmwfWindSpeed','owiWindSpeed_IPF', 'owiPreProcessing/estimWindDir200', 'owiPreProcessing/MixWindDir200', 'owiPreProcessing/CoefMixWindDir200', 'owiEcmwfWindDirection','owiWindDirection_mod1','owiWindSpeed_mod1','owiWindDirection_mod2','owiWindSpeed_mod2','owiWindDirection_mod3','owiWindSpeed_mod3','owiPreProcessing/FilterBinary200','owiLandFlag','owiIncidenceAngle')
        try:
            nc=netCDF4.Dataset(fname)
        except  Exception as e:
30
            logger.warning('skipping unreadable %s : %s' % (fname, str(e)))
ARCHER's avatar
ARCHER committed
31
32
33
34
35
36
37
38
39
40
            return None
        
        # read owiLon to get valid shape
        data['owiLon']=nc.variables['owiLon'][::]
        refDims=nc.variables['owiLon'].dimensions
        # read groups
        for groupName in [ None ] + list(nc.groups.keys()):
            if groupName:
                group=nc.groups[groupName]
                groupNamePref="%s/" % groupName
CEVAER's avatar
CEVAER committed
41
            else:
ARCHER's avatar
ARCHER committed
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
                group=nc
                groupNamePref=""
                
            groupAttrs=group.ncattrs()
            for varName in group.variables:
                var=group[varName]
                dims=var.dimensions
                if dims == refDims:
                    data["%s%s" % (groupNamePref,varName)] = var[::]
                else:
                    # some vars are 3D: 
                    try:
                        if dims[1:] == refDims:
                            values3D=var[::]
                            for islice in range(np.shape(values3D)[0]):
                                value=values3D[islice,:,:]
                                try:
                                    # try to find sliceName from varName and ncattrs
                                    sliceName=group.getncattr([ groupAttr for groupAttr in groupAttrs if groupAttr in varName ][0]).split(', ')[islice]
                                except:
                                    sliceName="%d" % islice
                                data["%s%s/%s" % (groupNamePref,varName, sliceName)]=value 
                    except:
                        pass
        
        nc.close
CEVAER's avatar
CEVAER committed
68
69
70
71
72
73
    
    return data

def readMeta(fname):
    data={}
    #list_var=('owiLon','owiLat','owiHeading','owiEcmwfWindSpeed','owiWindSpeed_IPF', 'owiPreProcessing/estimWindDir200', 'owiPreProcessing/MixWindDir200', 'owiPreProcessing/CoefMixWindDir200', 'owiEcmwfWindDirection','owiWindDirection_mod1','owiWindSpeed_mod1','owiWindDirection_mod2','owiWindSpeed_mod2','owiWindDirection_mod3','owiWindSpeed_mod3','owiPreProcessing/FilterBinary200','owiLandFlag','owiIncidenceAngle')
74
    attrs = ['long_name', 'units', 'valid_range', 'flag_values','flag_meanings', '_FillValue']
CEVAER's avatar
CEVAER committed
75
76
    try:
        nc=netCDF4.Dataset(fname)
77
    except  Exception as e:
78
        logger.warning('skipping unreadable %s : %s' % (fname, str(e)))
CEVAER's avatar
CEVAER committed
79
80
81
82
83
84
85
86
87
        return None
    
    # read globals attr
    for attr in nc.ncattrs():
        data[attr]=nc.getncattr(attr)
    
    refDims=nc.variables['owiLon'].dimensions
            
    # read groups
88
    for groupName in [ None ] + list(nc.groups.keys()):
CEVAER's avatar
CEVAER committed
89
90
91
92
93
94
95
96
97
98
99
        if groupName:
            group=nc.groups[groupName]
            groupNamePref="%s/" % groupName
        else:
            group=nc
            groupNamePref=""
            
        groupAttrs=group.ncattrs()
        for varName in group.variables:
            var=group[varName]
            data["%s%s" % (groupNamePref,varName)]={}
ARCHER's avatar
ARCHER committed
100
101
102
103
104
105
106
107
            for _attr in attrs:
                try:
                    data["%s%s" % (groupNamePref,varName)][_attr]=getattr(var,_attr)
                except:
                    data["%s%s" % (groupNamePref,varName)][_attr]=None
            # compatibility old version
            data["%s%s" % (groupNamePref,varName)]['description'] = data["%s%s" % (groupNamePref,varName)]['long_name']

CEVAER's avatar
CEVAER committed
108
109
110
111
112
113
114
115
116
117
118
119
120
            dims=var.dimensions
            if dims != refDims:
                
                # some vars are 3D: 
                try:
                    if dims[1:] == refDims:
                        for islice in range(var.shape[0]):
                            try:
                                # try to find sliceName from varName and ncattrs
                                sliceName=group.getncattr([ groupAttr for groupAttr in groupAttrs if groupAttr in varName ][0]).split(', ')[islice]
                            except:
                                sliceName="%d" % islice
                            data["%s%s/%s" % (groupNamePref,varName, sliceName)]={}
ARCHER's avatar
ARCHER committed
121
122
123
124
125
126
127
128
                            for _attr in attrs:
                                try:
                                    data["%s%s" % (groupNamePref,varName)][_attr]=getattr(var,_attr)
                                except:
                                    data["%s%s" % (groupNamePref,varName)][_attr]=None
                            # compatibility old version
                            data["%s%s" % (groupNamePref,varName)]['description'] = data["%s%s" % (groupNamePref,varName)]['long_name']

CEVAER's avatar
CEVAER committed
129
130
131
132
                except:
                    pass
    
    nc.close
ARCHER's avatar
ARCHER committed
133

CEVAER's avatar
CEVAER committed
134
135
    return data

136
137
138
139
def writeLight(fname,data,date,meta=None):
    """lightweigth file format writer    
    specs from https://github.com/Unidata/EC-netCDF-CF/blob/master/swath/swath.adoc#image-swath
    xarray reader exemple:
ARCHER's avatar
ARCHER committed
140
      * only root group:
141
        xarray.open_dataset(fname).transpose().squeeze() ####.rename_dims({'X':'x','Y':'y'})
ARCHER's avatar
ARCHER committed
142
143
144
145
      * all groups:
        xarray.merge([xarray.open_dataset(args.outfile),xarray.open_dataset(args.outfile,group='preprocessing')])
      
      
146
147
148
149
150
    """
    names_mapping = OrderedDict([
        ('owiLon'    , 'lon'),
        ('owiLat'    , 'lat'),
        ('owiWindSpeed' , 'wind_speed'),
151
        ('owiWindDirection' , 'wind_from_direction'),
ARCHER's avatar
ARCHER committed
152
153
        ('owiPreProcessing/wind_streaks_orientation' , 'wind_streaks_orientation'), 
        ('owiPreProcessing/wind_streaks_orientation_stddev' , 'wind_streaks_orientation_stddev'),
ARCHER's avatar
ARCHER committed
154
#        ('owiWindQuality' , 'quality_flag'),
155
        ('owiMask' , 'mask_flag'),
ARCHER's avatar
ARCHER committed
156
        ('owiWindFilter' , 'heterogeneity_mask'), 
157
158
159
160
161
162
        ('owiPreProcessing/FilterBinary_1' , 'heterogeneity_co_h'),
        ('owiPreProcessing/FilterBinary_2' , 'heterogeneity_co_m'),
        ('owiPreProcessing/FilterBinary_3' , 'heterogeneity_co_l'),
        ('owiPreProcessing/FilterBinary_1_cross' , 'heterogeneity_cr_h'),
        ('owiPreProcessing/FilterBinary_2_cross' , 'heterogeneity_cr_m'),
        ('owiPreProcessing/FilterBinary_3_cross' , 'heterogeneity_cr_l'),
ARCHER's avatar
ARCHER committed
163
164
        ('owiPreProcessing/ND_2','nrcs_detrend_co'), 
        ('owiPreProcessing/ND_2_cross' , 'nrcs_detrend_cross'),
ARCHER's avatar
ARCHER committed
165
        ('owiNrcs' , 'nrcs_co'),
ARCHER's avatar
ARCHER committed
166
        ('owiNrcs_cross' , 'nrcs_cross'),
ARCHER's avatar
ARCHER committed
167
        ('owiIncidenceAngle' , 'incidence_angle'),
ARCHER's avatar
ARCHER committed
168
        ('owiElevationAngle' , 'elevation_angle'),
169
        ])
ARCHER's avatar
ARCHER committed
170
171
172
173
174
175
176
177
178
    # global attrs (None mean no name translation )
    global_attrs = OrderedDict([
        ('sourceProduct' , None ), 
        ('missionName', None ), 
        ('polarisation', 'polarization' ), 
        ('footprint' , None ),
        ('l2ProcessingUtcTime' , None ),
        ('owiAlgorithmVersion', 'version' )
        ])
ARCHER's avatar
ARCHER committed
179
    attrs = ['long_name','units','valid_range','flag_values','flag_meanings']
180
181
182
    # overwrite meta
    meta_over = {
        'wind_speed' : {
ARCHER's avatar
ARCHER committed
183
184
            'units' : "m/s",
            'long_name' : "Ocean 10m Wind speed from co- and cross- polarization"
185
186
187
188
189
190
191
192
        },
        'lon' : {
            'standard_name' : 'longitude',
            'units'         : 'degrees_east'
        },
        'lat' : {
            'standard_name' : 'latitude',
            'units'         : 'degrees_north'
ARCHER's avatar
ARCHER committed
193
        },
194
195
        'wind_from_direction': {
            'long_name'     : 'Wind to direction (meteorologic convention)',
ARCHER's avatar
ARCHER committed
196
197
            'units'         : 'degrees',
#            'usage_level'   : "basic"
ARCHER's avatar
ARCHER committed
198
            }
199
    }
200
201
    # missing input will be forced to NaN/not implemented
    meta_not_implemented = {
202
        #'wind_streaks_orientation_stddev' : None
203
    }
204
205
206
207
208
    
    # translate dtype is only relevant for NETCDF4_CLASSIC
    translate_dtype = {
        np.dtype(np.uint8) : np.dtype(np.int16)
        }
209
    x,y = data['owiLon'].T.shape
210
    nc=netCDF4.Dataset(fname, 'w', format='NETCDF4_CLASSIC')
211
    nc.createDimension('time',1)
212
213
    nc.createDimension('y',y)
    nc.createDimension('x',x)
214
    nc.Conventions = 'CF-1.6'
ARCHER's avatar
ARCHER committed
215
    nc.title = 'SAR ocean surface wind field'
ARCHER's avatar
ARCHER committed
216
    nc.institution = 'IFREMER/CLS'
ARCHER's avatar
ARCHER committed
217
    nc.reference = 'Mouche Alexis, Chapron Bertrand, Knaff John, Zhao Yuan, Zhang Biao, Combot Clement (2019). Copolarized and Cross‐Polarized SAR Measurements for High‐Resolution Description of Major Hurricane Wind Structures: Application to Irma Category 5 Hurricane. Journal Of Geophysical Research-oceans, 124(6), 3905-3922. https://doi.org/10.1029/2019JC015056'
218
    nc.measurementDate = date.strftime("%Y-%m-%dT%H:%M:%SZ")
ARCHER's avatar
ARCHER committed
219
220
221
222
    for global_attr in global_attrs.keys():
        map_global_attr = global_attrs[global_attr]
        if map_global_attr is None:
            map_global_attr = global_attr
ARCHER's avatar
ARCHER committed
223
        try:
ARCHER's avatar
ARCHER committed
224
            setattr(nc, map_global_attr, meta[global_attr])
ARCHER's avatar
ARCHER committed
225
226
        except:
            logger.warning('skipping missing attr %s' % global_attr)
227
    time = nc.createVariable('time',np.int32,dimensions=('time')) # np.int32 induce 2038 bug, but is compatible with nc3
228
229
230
231
232
233
234
    time.standard_name = 'time'
    epoch = datetime.datetime(1970, 1, 1)
    time.units = "seconds since %s" % epoch.strftime("%Y-%m-%d %H:%M:%S")
    time.calendar = "gregorian"
    #time[::]  = int(date.timestamp())
    time[::] = (date - epoch).total_seconds() # works for python < 3.3
    for in_name in names_mapping.keys():
235
        out_name = names_mapping[in_name]
236
        if in_name in data:
ARCHER's avatar
ARCHER committed
237
            value = data[in_name].T
ARCHER's avatar
ARCHER committed
238
            fill_value = data[in_name].fill_value
ARCHER's avatar
ARCHER committed
239
            # transpose dims for correct netcdf ordering
240
241
242
243
            dtype = value.dtype
            if dtype in translate_dtype:
                dtype = translate_dtype[dtype]
            ncvar = nc.createVariable(out_name,dtype,dimensions=('time','y','x'),fill_value=fill_value)
244
        
ARCHER's avatar
ARCHER committed
245
            for _attr in attrs:           
246
                try:
ARCHER's avatar
ARCHER committed
247
248
249
250
251
                    try:
                        setattr(ncvar,_attr,meta_over[in_name][_attr])
                    except:
                        if meta[in_name][_attr] is not None:
                            setattr(ncvar,_attr, meta[in_name][_attr])
252
                except:
ARCHER's avatar
ARCHER committed
253
                    pass
254
            ncvar[::]=value.T
255
        else:
256
257
258
259
260
261
262
263
            if out_name in meta_not_implemented:
                logger.warning('Forcing not implemented var %s' % out_name)
                ncvar = nc.createVariable(out_name,np.float32,dimensions=('time','y','x'),fill_value=np.nan)
                setattr(ncvar,"long_name","Not yet impemented")
            
            else:
                logger.warning('skipping missing var %s' % in_name)
                continue
264
265
        if in_name not in ['owiLon','owiLat']:
            ncvar.coordinates="time lat lon"
ARCHER's avatar
ARCHER committed
266
        
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
        if out_name in meta_over:
            for tag in meta_over[out_name]:
                ncvar.setncattr(tag, meta_over[out_name][tag])
        
    #unable to set crs for gdal     
    #crs = nc.createVariable("crs",int)
    #crs.grid_mapping_name = "latitude_longitude"
    #crs.semi_major_axis = 6378137 
    #crs.inverse_flattening = 298.2572235604902
    #crs.longitude_of_prime_meridian = 0.0 
    
    
    
    nc.close()
    
    return 

ARCHER's avatar
ARCHER committed
284

285
286
287
288
289
def nclight2gridded(filein,fileout,crs_out=None):
    """
    if crs_out is None, each fileout will have his specific aeqd crs.
    if crs_out is not None, all fileout will share the same crs.
    """
ARCHER's avatar
ARCHER committed
290
291
292
293
294
295
296
297
    import numpy as np
    import xarray as xr
    import rioxarray as rxr
    import geo_shapely
    from pyproj.transformer import Transformer
    from shapely.geometry import Point
    from scipy import interpolate
    from scipy import ndimage
298
    import pyproj
ARCHER's avatar
ARCHER committed
299
    
300
301
302
303
304
305
306
307
308
309
310
311
312
    meta_over = {
        'x' : {
            'standard_name' : "projection_x_coordinate",
            'long_name' : 'Easting',
            'units'     : 'm'
            },
        'y' : {
            'standard_name' : "projection_y_coordinate",
            'long_name' : 'Northing',
            'units'     : 'm'
            },
        }
    
ARCHER's avatar
ARCHER committed
313
    ds_sw = xr.open_dataset(filein,mask_and_scale=False)
314
315
    x_center = ds_sw.x.size // 2
    y_center = ds_sw.y.size // 2
ARCHER's avatar
ARCHER committed
316
317
318
319
    
    lon_center = float(ds_sw.lon[0,y_center,x_center])
    lat_center = float(ds_sw.lat[0,y_center,x_center])

320
    # projection on aeqd . (will be intermediate crs if crs_out is not None)
ARCHER's avatar
ARCHER committed
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
    # get aeqd_crs at image center
    crs= geo_shapely.get_aeqd_crs(Point(lon_center,lat_center))
     
    # convert lon/lat in crs
    fromlonlat=Transformer.from_crs("epsg:4326",crs,always_xy=True) # always_xy=False is by default in pyproj 6+
    x_proj,y_proj=fromlonlat.transform(ds_sw.lon.values,ds_sw.lat.values)
    y_proj = y_proj.astype(float)
    x_proj = x_proj.astype(float)
    
    resolution = 1000 # in meters

    # compute grid size     
    x_size = int(np.ceil(np.max([-x_proj.min(),x_proj.max()])/resolution))*2+1
    y_size = int(np.ceil(np.max([-y_proj.min(),y_proj.max()])/resolution))*2+1
    
    
    # grid coordinates
    # do not use linspace : do it manualy, so computing inverse is easy
    #x = np.linspace(-(x_size-1)/2.,(x_size-1)/2.,num=x_size) * resolution
    #y = np.linspace(-(y_size-1)/2.,(y_size-1)/2.,num=y_size) * resolution
    x = (np.arange(x_size)-int((x_size-1)/2)) * resolution
    y = (np.arange(y_size)-int((y_size-1)/2)) * resolution
    # so given an x (resp y), in projection meters, we can have grid index with x / resolution
    x_idx = np.round((x_proj / resolution) + (x_size-1)/2 ).astype(int)
    y_idx = np.round((y_proj / resolution) + (y_size-1)/2 ).astype(int)
    
    # add x_proj and y_proj to original ds_sw, in valid coordinates (ie rounded at resolution)
    #ds_sw['x_proj'] = xr.DataArray(x_proj, dims=['time','Y', 'X'])
    #ds_sw['y_proj'] = xr.DataArray(y_proj, dims=['time','Y', 'X'])
    
    # get lon/lat for each grid cell
    #x_grid,y_grid = np.meshgrid(x,y)
    x_grid,y_grid = np.meshgrid(x,y)
    tolonlat=Transformer.from_crs(crs,"epsg:4326",always_xy=True)
    lon,lat=tolonlat.transform(x_grid,y_grid)
    
    # add time dim
    lon = lon.reshape([1] + list(lon.shape)).astype(np.float32)
    lat = lat.reshape([1] + list(lat.shape)).astype(np.float32)
    
    # generate dataarray
362
363
    da_lon = xr.DataArray(lon,coords=[ds_sw.time, y, x], dims=['time','y', 'x'])
    da_lat = xr.DataArray(lat,coords=[ds_sw.time, y, x], dims=['time','y', 'x'])
ARCHER's avatar
ARCHER committed
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
    
    # xarray dataset
    ds_gd = xr.Dataset({'lon' : da_lon,'lat' : da_lat})
    ds_gd.attrs = ds_sw.attrs
    
    # if resolution is too small, some gaps will occur in result
    # so we find gaps to be able to fill them later
    np_var = np.zeros(da_lon.shape[1:])
    # fill with 1 where the data will go
    np_var[y_idx,x_idx] = 1
    # mask is boolean grid where we need to interp
    mask = np_var == 0
    # xx and yy are all index
    xx, yy = np.meshgrid(np.arange(np_var.shape[1]), np.arange(np_var.shape[0]))
    # x_valid and y_valid are index of valid data
    #x_valid = xx[~mask]
    #y_valid = yy[~mask]
    # binary_closing and data to find isolated points that need interpolation
    mask_interp = ndimage.binary_closing(~mask) & mask
    # coordinates where interpolation is needed
    x_interp = xx[mask_interp]
    y_interp = yy[mask_interp]
    
    # variable loop
388
    for var in list((set(ds_sw.coords.keys()) | set(ds_sw.keys())).difference(set(['time']))):
ARCHER's avatar
ARCHER committed
389
390
        # set up dataarray. dtype is forced to np.float64, because we use Nan for interpolation
        #ds_sw[var].values.dtype
391
392
393
394
395
        try:
            fillvalue=ds_sw[var]._FillValue
        except:
            fillvalue=netCDF4.default_fillvals['f8']
        da_var = xr.DataArray(np.full(da_lon.shape, fillvalue, dtype=np.float64),coords=[ds_sw.time,y, x], dims=['time','y', 'x'])
ARCHER's avatar
ARCHER committed
396
397
398
399
        # fill dataarray
        da_var.values[0,y_idx,x_idx] = ds_sw[var].values
        # save nan, as they  are going to be interpolated
        nan_mask = np.isnan(da_var.values)
400
        da_var.values[nan_mask] = fillvalue # not really needed
ARCHER's avatar
ARCHER committed
401
402
        try:
            da_var.values[0,y_interp,x_interp] = np.nan
403
404
            da_var = da_var.ffill(dim='x',limit=2)
            da_var = da_var.ffill(dim='y',limit=2)
ARCHER's avatar
ARCHER committed
405
406
407
408
409
410
411
        except Exception as e:
            logging.warning("unable to fill nan on %s: %s" % (var,str(e)))
        # restore nan
        da_var.values[nan_mask] = np.nan
        # set original dtype
        da_var = da_var.astype(ds_sw[var].values.dtype)
        da_var.attrs = ds_sw[var].attrs
412
        da_var.rio.set_spatial_dims('x','y',inplace=True)
ARCHER's avatar
ARCHER committed
413
        ds_gd[var] = da_var
414
415

    encoding={}
416
    for var in list(set(ds_gd.coords.keys()) | set(ds_gd.keys())):
417
        encoding[var] = {}
418
419
        if var in meta_over:
            # python 3 only (merge dict )
ARCHER's avatar
ARCHER committed
420
421
422
            # ds_gd[var].attrs = { **ds_sw[var].attrs,**meta_over[var] }
            # python 2 compat:
            ds_gd[var].attrs.update(meta_over[var])
423
424
425
        if '_FillValue' not in ds_gd[var].attrs:
            # disable automatic fillvalue in fileout
            encoding[var].update({'_FillValue': None})
426
    
ARCHER's avatar
ARCHER committed
427
                
428
    ds_gd.rio.set_spatial_dims('x','y', inplace=True)
ARCHER's avatar
ARCHER committed
429
    ds_gd.rio.write_crs(crs, inplace=True)
430
    
431
432
    #encoding={'lat': {'_FillValue': None}, 'lon': {'_FillValue': None}}

433
434
435
436
    # final reprojection on crs_out
    if crs_out is not None:
        logger.info("reprojecting to final proj %s" % str(crs_out))
        ds_gd = ds_gd.rio.reproject(crs_out)
437
        if pyproj.CRS(crs_out).is_geographic :
438
439
440
441
442
443
444
445
446
447
448
449
            # crs out is lon/lat : rename x,y to lon,lat
            ds_gd = ds_gd.drop(['lon','lat']).rename({'x':'lon','y':'lat'})
            ds_gd['lon'].attrs = {
                "long_name" : "longitude",
                "standard_name" : "longitude",
                "units" : "degrees_east",
                }
            ds_gd['lat'].attrs = {
                "long_name" : "latitude",
                "standard_name" : "latitude",
                "units" : "degrees_north"
                }
450
451
            del encoding['x']
            del encoding['y']
452
            ds_gd.rio.set_spatial_dims('lon','lat', inplace=True)
453
            ds_gd.rio.write_crs(crs_out, inplace=True)
454
455
    
    ds_gd.to_netcdf(fileout,encoding=encoding,format='NETCDF4_CLASSIC')
456

CEVAER's avatar
CEVAER committed
457
458
459
def getAttribute(fname, attr):
    try:
        nc=netCDF4.Dataset(fname)
460
    except  Exception as e:
461
        logger.warning('skipping unreadable %s : %s' % (fname, str(e)))
CEVAER's avatar
CEVAER committed
462
463
464
465
466
467
468
469
470
        return None
    
    if attr in nc.ncattrs():
        return getattr(nc, attr)
    return None

def getFootprint(fname):
    try:
        nc=netCDF4.Dataset(fname)
471
    except  Exception as e:
472
        logger.warning('skipping unreadable %s : %s' % (fname, str(e)))
CEVAER's avatar
CEVAER committed
473
474
475
476
477
478
479
480
481
482
483
484
        return None
    
    if "footprint" in nc.ncattrs():
        return nc.footprint
    else:
        return None    
    
    
def getDateFromFname(fname):
    tmp=os.path.basename(fname).split('-')
    startDate=datetime.datetime.strptime(tmp[4],"%Y%m%dt%H%M%S")
    stopDate=datetime.datetime.strptime(tmp[5],"%Y%m%dt%H%M%S")
485
    date=startDate+old_div((stopDate-startDate),2)
CEVAER's avatar
CEVAER committed
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
    return date

def getDataAtCoord(data,lon,lat):
    import geo
    dataAtCoord={}
    #logger.info("lon : %s" % str(lon))
    landFlag=data['owiLandFlag']
    landFlag[landFlag>0]=1
    Mask = landFlag
    # Mask[Mask>0]=1
    # Mask1=data['owiLandFlag']
    # Mask2=np.logical_not(data['FilterBinary200'])
    lons = ma.array(data['owiLon'],mask=Mask)
    lats = ma.array(data['owiLat'],mask=Mask)
    # lons = ma.array(data['owiLon'],mask=Mask2)
    # lats = ma.array(data['owiLat'],mask=Mask2)
    # lons = data['owiLon']
    # lats = data['owiLat']
    dist,bearing=geo.haversine(lon, lat,lons, lats)
    mindist=np.min(dist)
    imindist=np.argmin(dist)
    if mindist > 20: #
508
        logger.warning("Buoy to far from sar image (%s km)" % mindist)
CEVAER's avatar
CEVAER committed
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
        return None
    data['dist_km']=mindist

    for key in data:
        try:
            dataAtCoord[key]=np.ravel(data[key])[imindist]
        except:
            dataAtCoord[key]=data[key]

        # convert dir in image conv to oceano convention
        # if ("Dir" in key) and (not key.startswith("DirectionModifie")) and (not key.startswith("Coef")) :
        if ("Dir200" in key) and (not key.startswith("DirectionModifie")) and (not key.startswith("Coef")) :
            if dataAtCoord[key] >180:
                dataAtCoord[key]=dataAtCoord[key]-180
            else :
                dataAtCoord[key]=dataAtCoord[key]+180
525
    return dataAtCoord
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561


# Convert L2C file system path into L2M filesystem path
# INPUT:
#     path : iterable or non iterable of a full L2C path (list of string or just string depending on input)
#     resolution : resolution of the L2M path needed
# Return:
#    A path or a list of path (same type as path parameter) containing the L2M paths
def L2CtoL2M(path, resolution):
    if type(path) == str:
        pathList = [path]
    else:
        pathList = path

    pathList = [p.replace("L2C", "L2M").replace("_CC_", "_CM_").replace("-cc-", "-cm-") for p in pathList]
    pathList = [p.split("/") for p in pathList]
    for p in pathList:
        filenameSplit = p[-1].split("-")
        # Set the field which indicates the file resolution to the needed resolution.
        # Adding leading zeros to keep field length
        filenameSplit[-2] = str(resolution).zfill(len(filenameSplit[-2]))
        p[-1] = "-".join(filenameSplit)
    pathList = ["/".join(p) for p in pathList]

    if type(path) == str:
        return pathList[0]
    else:
        return pathList


# Convert L2C file system path into L2C/L2M light filesystem path
# INPUT:
#     path (list of string or string) : iterable or non iterable of a full L2C path
#     resolution (int) : resolution of the L2M path needed. If not provided, will return the L2CLight
# Return:
#    A path or a list of path (same type as path parameter) containing the L2CLight/L2MLight paths
ARCHER's avatar
ARCHER committed
562
def L2CtoLight(path, resolution=None , mode='sw'):
563
564
565
566
567
568
569
570
571
572
573
574
575
    if type(path) == str:
        pathList = [path]
    else:
        pathList = path

    pathsOut = []
    for p in pathList:
        base = os.path.dirname(p)
        if resolution:
            addPath = "post_processing/nclight_L2M"
        else:
            addPath = "post_processing/nclight"
        newPath = os.path.join(base, addPath)
ARCHER's avatar
ARCHER committed
576
        newNcFile = os.path.splitext(os.path.basename(p))[0] + "_%s.nc" % mode
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
        if resolution:
            newNcFile = newNcFile.replace("-cc-", "-cm-")
            filenameSplit = newNcFile.split("-")
            # Set the field which indicates the file resolution to the needed resolution.
            # Adding leading zeros to keep field length
            filenameSplit[-2] = str(resolution).zfill(len(filenameSplit[-2]))
            newNcFile = "-".join(filenameSplit)
        fullNcPath = os.path.join(newPath, newNcFile)
        pathsOut.append(fullNcPath)

    if type(path) == str:
        return pathsOut[0]
    else:
        return pathsOut