Commit f2d7ba78 authored by ARCHER's avatar ARCHER

nclight gridded can share same global grid

parent 7b22b8ad
......@@ -23,7 +23,9 @@ if __name__ == "__main__":
parser.add_argument('outfile')
args = parser.parse_args()
logger.info('%s -> %s' % (args.infile,args.outfile))
owi.nclight2gridded(args.infile,args.outfile)
crs_out = "epsg:4326"
#crs_out = None
owi.nclight2gridded(args.infile,args.outfile,crs_out= crs_out)
# 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
......@@ -287,7 +287,11 @@ def writeLight(fname,data,date,meta=None):
return
def nclight2gridded(filein,fileout):
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.
"""
import numpy as np
import xarray as xr
import rioxarray as rxr
......@@ -296,6 +300,7 @@ def nclight2gridded(filein,fileout):
from shapely.geometry import Point
from scipy import interpolate
from scipy import ndimage
import pyproj
meta_over = {
'x' : {
......@@ -317,6 +322,7 @@ def nclight2gridded(filein,fileout):
lon_center = float(ds_sw.lon[0,y_center,x_center])
lat_center = float(ds_sw.lat[0,y_center,x_center])
# projection on aeqd . (will be intermediate crs if crs_out is not None)
# get aeqd_crs at image center
crs= geo_shapely.get_aeqd_crs(Point(lon_center,lat_center))
......@@ -417,7 +423,31 @@ def nclight2gridded(filein,fileout):
ds_gd.rio.set_spatial_dims('x','y', inplace=True)
ds_gd.rio.write_crs(crs, inplace=True)
ds_gd.to_netcdf(fileout,format='NETCDF4_CLASSIC')
encoding=None
# 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)
if pyproj.CRS(crs_out).is_geographic :
# 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"
}
encoding = {'lat': {'_FillValue': None}, 'lon': {'_FillValue': None}}
ds_gd.rio.set_spatial_dims('lon','lat', inplace=True)
ds_gd.rio.write_crs(crs, inplace=True)
ds_gd.to_netcdf(fileout,encoding=encoding,format='NETCDF4_CLASSIC')
def getAttribute(fname, attr):
try:
......
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