Commit 89c6a9f2 authored by ARCHER's avatar ARCHER 💬

gridded nclight

parent 84d8392b
#!/usr/bin/env python
import sys
import owi
import argparse
import sys
import netCDF4
import datetime
import logging
logging.basicConfig()
logger = logging.getLogger(__name__)
if sys.gettrace():
logger.setLevel(logging.DEBUG)
else:
logger.setLevel(logging.INFO)
if __name__ == "__main__":
description = """Converts nclight file to gridded nc format"""
parser = argparse.ArgumentParser(description = description)
parser.add_argument('infile')
parser.add_argument('outfile')
args = parser.parse_args()
logger.info('%s -> %s' % (args.infile,args.outfile))
owi.nclight2gridded(args.infile,args.outfile)
# gdaldem color-relief NETCDF:"/tmp/test_gd.nc":wind_speed wind_speed.txt /tmp/testowi1.tif -alpha
# gdal_translate -of GTiff NETCDF:"/tmp/test_gd.nc":wind_speed /tmp/testowi.tif
......@@ -259,10 +259,118 @@ def writeLight(fname,data,date,meta=None):
return
#def readLight(fname):
# import xarray
# ds = xarray.open_dataset(fname)
# return ds
def nclight2gridded(filein,fileout):
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
ds_sw = xr.open_dataset(filein,mask_and_scale=False)
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])
# 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
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})
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
for var in ds_sw.keys():
# 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'])
# 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)
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
ds_gd[var] = da_var
ds_gd.rio.set_spatial_dims('X','Y', inplace=True)
ds_gd.rio.write_crs(crs, inplace=True)
ds_gd.to_netcdf(fileout)
def getAttribute(fname, attr):
try:
......@@ -369,7 +477,7 @@ def L2CtoL2M(path, resolution):
# 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
def L2CtoLight(path, resolution=None):
def L2CtoLight(path, resolution=None , mode='sw'):
if type(path) == str:
pathList = [path]
else:
......@@ -383,7 +491,7 @@ def L2CtoLight(path, resolution=None):
else:
addPath = "post_processing/nclight"
newPath = os.path.join(base, addPath)
newNcFile = os.path.splitext(os.path.basename(p))[0] + "_sw.nc"
newNcFile = os.path.splitext(os.path.basename(p))[0] + "_%s.nc" % mode
if resolution:
newNcFile = newNcFile.replace("-cc-", "-cm-")
filenameSplit = newNcFile.split("-")
......
......@@ -10,7 +10,7 @@ setup(name='owi',
setup_requires=['setuptools_scm'],
packages=['owi'],
zip_safe=False,
scripts=['bin/owi2nclight.py'],
scripts=['bin/owi2nclight.py','bin/nclight2gridded.py'],
install_requires=[
'numpy', 'netCDF4'
])
......
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