Commit 7e1a5df6 authored by ARCHER's avatar ARCHER

added not implemented fields, renamed dims ('X','Y') to ('x','y')

parent 90d3bc0a
......@@ -138,7 +138,7 @@ def writeLight(fname,data,date,meta=None):
specs from https://github.com/Unidata/EC-netCDF-CF/blob/master/swath/swath.adoc#image-swath
xarray reader exemple:
* only root group:
xarray.open_dataset(fname).transpose().squeeze().rename_dims({'X':'x','Y':'y'})
xarray.open_dataset(fname).transpose().squeeze() ####.rename_dims({'X':'x','Y':'y'})
* all groups:
xarray.merge([xarray.open_dataset(args.outfile),xarray.open_dataset(args.outfile,group='preprocessing')])
......@@ -160,7 +160,9 @@ def writeLight(fname,data,date,meta=None):
# ('owiPreProcessing/FilterBinary_1_cross' , 'heterogeneity_cr_h'),
# ('owiPreProcessing/FilterBinary_2_cross' , 'heterogeneity_cr_m'),
# ('owiPreProcessing/FilterBinary_3_cross' , 'heterogeneity_cr_l'),
('owiPreProcessing/FilterBinary_2' , 'heterogeneity_mask'),
('owiPreProcessing/FilterBinary_2' , 'heterogeneity_mask'), # TODO: waiting for merged filter
('owiPreProcessing/wind_streaks_orientation_stddev' , 'wind_streaks_orientation_stddev'),
('owiPreProcessing/estimWindDir_2' , 'wind_streaks_orientation'), # TODO: should be owiPreProcessing/wind_streaks_orientation
('owiNrcs' , 'nrcs_co'),
('owiNrcs_cross' , 'nrcs_cr'),
('owiElevationAngle' , 'elevation_angle'),
......@@ -189,11 +191,15 @@ def writeLight(fname,data,date,meta=None):
# 'usage_level' : "basic"
}
}
X,Y = data['owiLon'].T.shape
# missing input will be forced to NaN/not implemented
meta_not_implemented = {
'wind_streaks_orientation_stddev' : None
}
x,y = data['owiLon'].T.shape
nc=netCDF4.Dataset(fname, 'w', format='NETCDF4')
nc.createDimension('time',1)
nc.createDimension('Y',Y)
nc.createDimension('X',X)
nc.createDimension('y',y)
nc.createDimension('x',x)
nc.Conventions = 'CF-1.6'
nc.title = 'SAR ocean surface wind field'
nc.institution = 'IFREMER/CLS'
......@@ -213,15 +219,15 @@ def writeLight(fname,data,date,meta=None):
#time[::] = int(date.timestamp())
time[::] = (date - epoch).total_seconds() # works for python < 3.3
for in_name in names_mapping.keys():
out_name = names_mapping[in_name]
if in_name in data:
value = data[in_name].T
out_name = names_mapping[in_name]
try:
fill_value=data[in_name].fill_value
except:
fill_value=None
# transpose dims for correct netcdf ordering
ncvar = nc.createVariable(out_name,value.dtype,dimensions=('time','Y','X'),fill_value=fill_value)
ncvar = nc.createVariable(out_name,value.dtype,dimensions=('time','y','x'),fill_value=fill_value)
for _attr in attrs:
try:
......@@ -232,18 +238,19 @@ def writeLight(fname,data,date,meta=None):
setattr(ncvar,_attr, meta[in_name][_attr])
except:
pass
ncvar[::]=value.T
else:
logger.warning('skipping missing var %s' % in_name)
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
if in_name not in ['owiLon','owiLat']:
ncvar.coordinates="time lat lon"
if in_name not in ['owiWindQuality']:
ncvar[::]=value.T
else:
# known bad vars are set to _FillValue
ncvar[::]=ncvar._FillValue
#ncvar.grid_mapping = "crs"
if out_name in meta_over:
for tag in meta_over[out_name]:
ncvar.setncattr(tag, meta_over[out_name][tag])
......@@ -272,9 +279,22 @@ def nclight2gridded(filein,fileout):
from scipy import interpolate
from scipy import ndimage
meta_over = {
'x' : {
'standard_name' : "projection_x_coordinate",
'long_name' : 'Easting',
'units' : 'm'
},
'y' : {
'standard_name' : "projection_y_coordinate",
'long_name' : 'Northing',
'units' : 'm'
},
}
ds_sw = xr.open_dataset(filein,mask_and_scale=False)
x_center = ds_sw.X.size // 2
y_center = ds_sw.Y.size // 2
x_center = ds_sw.x.size // 2
y_center = ds_sw.y.size // 2
lon_center = float(ds_sw.lon[0,y_center,x_center])
lat_center = float(ds_sw.lat[0,y_center,x_center])
......@@ -320,8 +340,8 @@ def nclight2gridded(filein,fileout):
lat = lat.reshape([1] + list(lat.shape)).astype(np.float32)
# generate dataarray
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'])
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'])
# xarray dataset
ds_gd = xr.Dataset({'lon' : da_lon,'lat' : da_lat})
......@@ -346,21 +366,19 @@ def nclight2gridded(filein,fileout):
y_interp = yy[mask_interp]
# variable loop
for var in ds_sw.keys():
for var in list((set(ds_sw.coords.keys()) | set(ds_sw.keys())).difference(set(['time']))):
# set up dataarray. dtype is forced to np.float64, because we use Nan for interpolation
#ds_sw[var].values.dtype
da_var = xr.DataArray(np.full(da_lon.shape, ds_sw[var]._FillValue, dtype=np.float64),coords=[ds_sw.time,y, x], dims=['time','Y', 'X'])
da_var = xr.DataArray(np.full(da_lon.shape, ds_sw[var]._FillValue, dtype=np.float64),coords=[ds_sw.time,y, x], dims=['time','y', 'x'])
# fill dataarray
# FIXME : ds_sw[var].values must not contain Nan : they will be interpolated
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)
da_var.values[nan_mask] = ds_sw[var]._FillValue # not really needed
try:
da_var.values[0,y_interp,x_interp] = np.nan
da_var = da_var.ffill(dim='X',limit=2)
da_var = da_var.ffill(dim='Y',limit=2)
da_var = da_var.ffill(dim='x',limit=2)
da_var = da_var.ffill(dim='y',limit=2)
except Exception as e:
logging.warning("unable to fill nan on %s: %s" % (var,str(e)))
# restore nan
......@@ -368,9 +386,17 @@ def nclight2gridded(filein,fileout):
# set original dtype
da_var = da_var.astype(ds_sw[var].values.dtype)
da_var.attrs = ds_sw[var].attrs
da_var.rio.set_spatial_dims('x','y',inplace=True)
ds_gd[var] = da_var
for var in list(set(ds_gd.coords.keys()) | set(ds_gd.keys())):
if var in meta_over:
# python 3 only (merge dict )
ds_gd[var].attrs = { **ds_sw[var].attrs,**meta_over[var] }
ds_gd.rio.set_spatial_dims('X','Y', inplace=True)
ds_gd.rio.set_spatial_dims('x','y', inplace=True)
ds_gd.rio.write_crs(crs, inplace=True)
ds_gd.to_netcdf(fileout)
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment