From c6cc0976cd446fc91fdac56de14b72da9b7bfcd6 Mon Sep 17 00:00:00 2001 From: BaptisteVandecrux Date: Thu, 8 Sep 2022 14:08:51 +0200 Subject: [PATCH] Small IO update to compare with pySICEv2.0 --- sice_f.py | 437 +++++++++++++++++++++++++--------------------------- sice_io.py | 47 +++--- sice_lib.py | 42 ++--- 3 files changed, 244 insertions(+), 282 deletions(-) diff --git a/sice_f.py b/sice_f.py index ea1324f..bdf02bc 100644 --- a/sice_f.py +++ b/sice_f.py @@ -2,13 +2,10 @@ """ Created on Fri Oct 25 11:35:10 2019 -This script reads a set of input tif files and outputs the csv needed to run sice.f +tip list: + %matplotlib inline + %matplotlib qt -You will need to update - InputFoder - the path to water_vod and ozone_vod - the output path - @author: bav@geus.dk """ @@ -17,245 +14,229 @@ import time import os import sys -import bav_lib as bl +#import bav_lib as bl import pandas as pd - +import rioxarray +import shutil +import glob +import subprocess +InputFolder = './data/2019-06-14/' +OutputFolder = './data/2019-06-14/fortran/' +ProcessingFolder = './fortran/' +# %% start_time = time.time() +InputFolder = InputFolder + "/" +Oa01 = rio.open(InputFolder + "r_TOA_01.tif") +meta = Oa01.meta -InputFolder = sys.argv[1] + '/' -# InputFolder = 'out/SICE_2020_py1.4/' - -print(os.getcwd()) - +def WriteOutput(var, var_name, in_folder): + # this functions write tif files based on a model file, here "Oa01" + # opens a file for writing + with rio.open(in_folder + var_name + ".tif", "w+", **meta) as dst: + dst.write(var.astype("float32"), 1) -#%% turning geotiffs into sice input -print("Reading input ") -print(sys.argv[1]) - -#%% input text file -if os.path.isfile(sys.argv[1]): - InputFolder = os.path.dirname(sys.argv[1])+'/' - # data_in = pd.read_csv('validation/data/S3_PROMICE.csv') - data_in = pd.read_csv(sys.argv[1]) - toa = np.expand_dims(data_in[[c for c in data_in.columns if c.find('reflec')>=0]].to_numpy().transpose(), axis=2) - - - ozone = np.expand_dims(data_in['total_ozone'], axis=1) - water = np.expand_dims(data_in['total_columnar_water_vapour'], axis=1) - sza = np.expand_dims(data_in['sza'], axis=1) - saa = np.expand_dims(data_in['saa'], axis=1) - vza = np.expand_dims(data_in['vza'], axis=1) - vaa = np.expand_dims(data_in['vaa'], axis=1) - height = np.expand_dims(data_in['altitude'], axis=1) - - sza[np.isnan(toa[0,:,:])] = np.nan - saa[np.isnan(toa[0,:,:])] = np.nan - vza[np.isnan(toa[0,:,:])] = np.nan - vaa[np.isnan(toa[0,:,:])] = np.nan - -#%% ========= input tif =============== -elif os.path.isdir(sys.argv[1]): - print('\n tiff input') - InputFolder = sys.argv[1] + '/' - Oa01 = rio.open(InputFolder+'r_TOA_01.tif') - meta = Oa01.meta - - def WriteOutput(var,var_name,in_folder): - # this functions write tif files based on a model file, here "Oa01" - # opens a file for writing - with rio.open(in_folder+var_name+'.tif', 'w+', **meta) as dst: - dst.write(var.astype('float32'),1) - - toa = np.tile(Oa01.read(1)*np.nan, (21,1,1)) - - for i in range(21): - dat = rio.open((InputFolder+'r_TOA_'+str(i+1).zfill(2)+'.tif')) - toa[i,:,:] = dat.read(1) - - ozone = rio.open(InputFolder+'O3.tif').read(1) - water = rio.open(InputFolder+'WV.tif').read(1) - sza = rio.open(InputFolder+'SZA.tif').read(1) - saa = rio.open(InputFolder+'SAA.tif').read(1) - vza = rio.open(InputFolder+'OZA.tif').read(1) - vaa = rio.open(InputFolder+'OAA.tif').read(1) - height = rio.open(InputFolder+'height.tif').read(1) - - sza[np.isnan(toa[0,:,:])] = np.nan - saa[np.isnan(toa[0,:,:])] = np.nan - vza[np.isnan(toa[0,:,:])] = np.nan - vaa[np.isnan(toa[0,:,:])] = np.nan - -# ns,alat,alon,sza,vza,saa,vaa,height, (toa(iks),iks=1,21),OZON,WATER -olci_toa = np.vstack((np.arange(1,len(sza.flatten())+1), - sza.flatten()*np.nan, - sza.flatten()*np.nan, - sza.flatten(), - vza.flatten(), - saa.flatten(), - vaa.flatten(), - height.flatten())) +toa = np.tile(Oa01.read(1) * np.nan, (21, 1, 1)) for i in range(21): - olci_toa = np.vstack((olci_toa, toa[i,:,:].flatten())) - + dat = rio.open((InputFolder + "r_TOA_" + str(i + 1).zfill(2) + ".tif")) + toa[i, :, :] = dat.read(1) + +ozone = rio.open(InputFolder + "O3.tif").read(1) +water = rio.open(InputFolder + "WV.tif").read(1) +sza = rio.open(InputFolder + "SZA.tif").read(1) +saa = rio.open(InputFolder + "SAA.tif").read(1) +vza = rio.open(InputFolder + "OZA.tif").read(1) +vaa = rio.open(InputFolder + "OAA.tif").read(1) +height = rio.open(InputFolder + "height.tif").read(1) + +sza[np.isnan(toa[0, :, :])] = np.nan +saa[np.isnan(toa[0, :, :])] = np.nan +vza[np.isnan(toa[0, :, :])] = np.nan +vaa[np.isnan(toa[0, :, :])] = np.nan + +olci_toa = np.vstack( + ( + np.arange(1, len(sza.flatten()) + 1), # pixel number_x + np.arange(1, len(sza.flatten()) + 1), # pixel number_y + np.arange(1, len(sza.flatten()) + 1), # latitude + np.arange(1, len(sza.flatten()) + 1), # longitude + sza.flatten(), # solar zenith angle + saa.flatten(), # soalr azimuthal angle + vza.flatten(), # viewing zenith angle + vaa.flatten(), # viewing azimuth angle + ) +) +for i in range(21): + olci_toa = np.vstack((olci_toa, toa[i, :, :].flatten())) +olci_toa = np.vstack((olci_toa, height.flatten())) olci_toa = np.vstack((olci_toa, ozone.flatten())) -olci_toa = np.vstack((olci_toa, water.flatten())) -olci_toa=olci_toa.T +olci_toa = olci_toa.T olci_toa_save = olci_toa -ind_good_pixels = np.logical_not(np.isnan(toa[0,:,:]).flatten()) -olci_toa = olci_toa[ind_good_pixels,:] +ind_good_pixels = np.logical_not(np.isnan(toa[0, :, :]).flatten()) +olci_toa = olci_toa[ind_good_pixels, :] olci_toa[np.isnan(olci_toa)] = 999 -OutputFolder= InputFolder+'fortran/' -#os.mkdir(OutputFolder) -print('\nInput file saved: '+OutputFolder+'olci_toa_newformat.dat') -np.savetxt(OutputFolder+'olci_toa_newformat.dat', X=olci_toa, delimiter='\t', \ - fmt='%i '+ '%10.5f '*6 + '%i '+ '%10.5f'*23) - -#%% You can now run sice.f -# print('bash ./fortran/sice_f.sh -i '+OutputFolder) -# os.system('bash ./fortran/sice_f.sh -i '+OutputFolder) - -cmnd = 'bash ./fortran/sice_f.sh -i '+OutputFolder -print(cmnd) -import subprocess - -try: - output = subprocess.check_output( - cmnd, stderr=subprocess.STDOUT, shell=True, universal_newlines=True) -except subprocess.CalledProcessError as exc: - print("Status : FAIL", exc.returncode, exc.output) -else: - print("Output: \n{}\n".format(output)) - -#%% Loading result from sice_debug.f -# OUTPUT files: -# file= 'spherical_albedo.dat' ) ns,ndate(3),alat,alon,(answer(i),i=1,21),isnow -# file= 'planar_albedo.dat' ) ns,ndate(3),alat,alon,(rp(i),i=1,21),isnow -# file= 'boar.dat' ) ns,ndate(3),alat,alon,(refl(i),i=1,21),isnow -# file= 'size.dat' ) ns,ndate(3),alat,alon,D,area,al,r0, andsi,andbi,indexs,indexi,indexd,isnow -# file= 'impurity.dat' ) ns,ndate(3),alat,alon,ntype,conc,bf,bm,thv, toa(1),isnow -# file= 'bba.dat' ) ns,ndate(3),alat,alon,rp3,rs3, isnow -# file= 'bba_alex_reduced.dat') ns,ndate(3),rp3,isnow -# file= 'notsnow.dat') ns,ndate(3),icloud,iice -# file='retrieved_O3.dat') ns,alat,alon,BXXX,totadu,deltak,sza,vza,amf - -if os.path.isfile(sys.argv[1]): - impurity = pd.read_csv(OutputFolder+'impurity.dat', - names=['ns','ndate','alat','alon','ntype','conc', - 'bf','bm','thv', 'toa(1)','isnow'], - sep=' ', skipinitialspace=True, header=None) - size = pd.read_csv(OutputFolder+'size.dat', - names=['ns','ndate','alat','alon','D','area','al', - 'r0', 'andsi','andbi','indexs','indexi','indexd','isnow'], - sep=' ', skipinitialspace=True, header=None) - bba = pd.read_csv(OutputFolder+'bba.dat', - names=['ns','ndate','alat','alon','rp3','rs3','isnow'], - sep=' ', skipinitialspace=True, header=None) - spherical_albedo = pd.read_csv(OutputFolder+'spherical_albedo.dat', - sep=' ', skipinitialspace=True, header=None).values - planar_albedo = pd.read_csv(OutputFolder+'planar_albedo.dat', - sep=' ', skipinitialspace=True, header=None).values - data_out = data_in - data_out['grain_diameter']=np.nan - data_out['snow_specific_area']=np.nan - data_out['al']=np.nan - data_out['r0']=np.nan - data_out['diagnostic_retrieval']=np.nan - data_out['conc']=np.nan - data_out['albedo_bb_planar_sw']=np.nan - data_out['albedo_bb_spherical_sw']=np.nan +# os.mkdir(OutputFolder) +print("\nInput file saved: " + OutputFolder + "input.dat") +np.savetxt( + OutputFolder + "input.dat", + X=olci_toa, + delimiter="\t", + fmt="%i "*2 + "%10.5f " * 29, +) - data_out['grain_diameter'][size.ns-1] = size.D - # data_out.loc[size.ns-1, ['grain_diameter']]=size.D - data_out.loc[size.ns-1,['snow_specific_area']]=size.area - data_out.loc[size.ns-1,['al']]=size.al - data_out.loc[size.ns-1,['r0']]=size.r0 - data_out.loc[size.ns-1,['diagnostic_retrieval']]=bba.isnow - data_out.loc[size.ns-1,['conc']]=impurity.conc - data_out['albedo_bb_planar_sw'][size.ns-1]=bba.rp3 - data_out['albedo_bb_spherical_sw']=np.nan - data_out['albedo_bb_spherical_sw'][size.ns-1]=bba.rs3 - - # data_out.iloc[size.ns-1, data_out.columns.get_loc('grain_diameter')]=size.D - # data_out.iloc[size.ns-1, data_out.columns.get_loc('snow_specific_area')]=size.area - # data_out.iloc[size.ns-1, data_out.columns.get_loc('al')]=size.al - # data_out.iloc[size.ns-1, data_out.columns.get_loc('r0')]=size.r0 - # data_out.iloc[size.ns-1, data_out.columns.get_loc('diagnostic_retrieval')]=bba.isnow - # data_out.iloc[size.ns-1, data_out.columns.get_loc('conc')]=impurity.conc - # data_out.iloc[size.ns-1, data_out.columns.get_loc('albedo_bb_planar_sw')]=bba.rp3 - # data_out.iloc[size.ns-1, data_out.columns.get_loc('albedo_bb_spherical_sw')]=bba.rs3 +#%% You can now run sice.f +shutil.copy(OutputFolder+'input.dat', ProcessingFolder[:-1]) +os.chdir('./fortran') +print('Running sice.f') +#subprocess.run('gfortran ./version\ 6.1\ 2022/s.f -o ./sice.exe', shell=True) +subprocess.check_call('./sice.exe') +for file in glob.glob(r'*.dat'): + if file == 'thv.dat': + continue + if file == 'input.dat': + continue - - for i in np.append(np.arange(11), np.arange(15,21)): - # for i in np.arange(21): - data_out['albedo_spectral_spherical_'+str(i+1).zfill(2)]=np.nan - data_out.loc[size.ns-1,['albedo_spectral_spherical_'+str(i+1).zfill(2)]]=spherical_albedo[:,4+i] + print('Moving '+file) + shutil.move(file, '../'+OutputFolder+'/'+file) +os.chdir('..') + + +#%% Converting text output to geotiff +Oa01 = rio.open(InputFolder + "r_TOA_01.tif") +meta = Oa01.meta +with rio.Env(): + meta.update(compress="DEFLATE") + +def output_sice_f(file_name, var_name, var_id): + sice_out = pd.read_csv(file_name, delim_whitespace=True, header=None).values + ind_orig = np.arange(1, len(Oa01.read(1).flatten()) + 1) + var = ind_orig * np.nan + ind_pix = sice_out[:, 1].astype(int) + var[ind_pix - 1] = sice_out[:, var_id] + var_mat = np.reshape(var, np.shape(Oa01.read(1))) + var_mat[var_mat == 999] = np.nan + with rio.open(OutputFolder + var_name + ".tif", "w+", **meta) as dst: + dst.write(var_mat.astype("float32"), 1) + +output_sice_f(OutputFolder + "snow_parameters.dat", "isnow", 3) +output_sice_f(OutputFolder + "snow_parameters.dat", "factor", 4) +output_sice_f(OutputFolder + "snow_parameters.dat", "grain_diameter", 5) +output_sice_f(OutputFolder + "snow_parameters.dat", "snow_specific_area", 6) +output_sice_f(OutputFolder + "snow_parameters.dat", "al", 7) +output_sice_f(OutputFolder + "snow_parameters.dat", "r0", 8) +output_sice_f(OutputFolder + "snow_parameters.dat", "bm", 10) +output_sice_f(OutputFolder + "snow_parameters.dat", "polut", 11) +output_sice_f(OutputFolder + "snow_parameters.dat", "albedo_bb_planar_sw", 16) +output_sice_f(OutputFolder + "snow_parameters.dat", "albedo_bb_spherical_sw", 19) + +i = 0 +output_sice_f( + OutputFolder + "spectral_spherical_albedo.dat", + "alb_sph_" + str(i + 1).zfill(2), + 5 + i, +) +output_sice_f( + OutputFolder + "spectral_spherical_albedo_solved.dat", + "alb_sph_" + str(i + 1).zfill(2) + '_solved', + 5 + i, +) +output_sice_f( + OutputFolder + "spectral_plane_albedo.dat", + "alb_pl_" + str(i + 1).zfill(2), + 5 + i, +) + +output_sice_f( + OutputFolder + "spectral_plane_albedo_solved.dat", + "alb_pl_" + str(i + 1).zfill(2)+ '_solved', + 5 + i, +) +output_sice_f( + OutputFolder + "spectral_bOAR.dat", + "BOAR_" + str(i + 1).zfill(2), + 5 + i, +) + +output_sice_f( + OutputFolder + "spectral_BOAR_solved.dat", + "bOAR_" + str(i + 1).zfill(2)+ '_solved', + 5 + i, +) + + +# for i in range(21): +# output_sice_f( +# OutputFolder + "spherical_albedo.dat", +# "albedo_spectral_spherical_" + str(i + 1).zfill(2), +# 4 + i, +# ) +# output_sice_f( +# OutputFolder + "planar_albedo.dat", +# "albedo_spectral_planar_" + str(i + 1).zfill(2), +# 4 + i, +# ) +# format snow_parameters.dat: +# "j", "alat", "alon", "NCLASS", +# "factor", "diam", "ssa", "dlina", "rv", "aload1", +# "powe", "polut", "eff","absor1","absef660","absor1000", +# "rsw", "rvis","rnir", "rsws", "rviss","rnirs", +# "andbi","andsi","ratka", "NPOLE", +# "NBARE", "NSNOW", "sza","vza","raa","toa(1)","toa(21)", +# "tocos","akozon","difka","cv1","cv2" +# format spectral_spherical_albedo.dat: +# WRITE(1001,*) j,alat,alon, NCLASS,factor,(albs(ir),ir=1,21) + +#%% plotting oulooks +import xarray as xr +import rioxarray +import matplotlib.pyplot as plt +#%matplotlib qt +import numpy as np +data_folder = 'data/2019-06-14/' +code_ver1 = 'pySICEv1.6' +code_ver2 = 'pySICEv2.0' +folder1 = data_folder+code_ver1+'/' +folder2 = data_folder+code_ver2+'/' +var_list = ['isnow', "grain_diameter", "snow_specific_area", + "r0", "albedo_bb_planar_sw", "albedo_bb_spherical_sw"] +# var_list = ["isnow", "alb_sph_01", "alb_sph_01_solved", "alb_pl_01", "alb_pl_01_solved", "BOAR_01", "BOAR_01_solved"] +plt.close('all') +for var in var_list: + ds_f = rioxarray.open_rasterio(folder1+var+'.tif').squeeze() + ds_p = rioxarray.open_rasterio(folder2+var+'.tif').squeeze().interp_like(ds_f, method ='nearest') + + if var == 'polut': + ds_f = np.log10(ds_f) + ds_p = np.log10(ds_p) - for i in np.append(np.arange(11), np.arange(15,21)): - data_out['rBRR_'+str(i+1).zfill(2)]=np.nan - data_out.loc[size.ns-1,['rBRR_'+str(i+1).zfill(2)]]=planar_albedo[:,4+i] + if var in ['albedo_bb_planar_sw', 'albedo_bb_spherical_sw']: + ds_f = ds_f.where(ds_f>0).where(ds_f<2) + if var == 'isnow': + ds_f = ds_f.where(ds_f>=0).where(ds_f<4) - data_out.to_csv(sys.argv[1][:-4]+'_fortran_out.csv') - print('\nOutput: '+ sys.argv[1][:-4] + '_fortran_out.csv') + + vmin = min(ds_f.min(), ds_p.min()) + vmax = max(ds_f.max(), ds_p.max()) + fig, ax = plt.subplots(2,2, figsize=(15,15)) + ax=ax.flatten() -# ========= input tif =============== -elif os.path.isdir(sys.argv[1]): - Oa01 = rio.open(InputFolder+'r_TOA_01.tif') - meta = Oa01.meta - with rio.Env(): - meta.update(compress='DEFLATE') - def output_sice_f(file_name,var_name,var_id): - sice_out = np.loadtxt(file_name) - ind_orig = np.arange(1,len(Oa01.read(1).flatten())+1) - var = ind_orig*np.nan - ind_pix = sice_out[:,0].astype(int) - var[ind_pix-1] = sice_out[:,var_id] - var_mat = np.reshape(var,np.shape(Oa01.read(1))) - var_mat[var_mat==999] = np.nan - with rio.open(OutputFolder+var_name+'.tif', 'w+', **meta) as dst: - dst.write(var_mat.astype('float32'),1) - output_sice_f(OutputFolder+"bba.dat",'albedo_bb_planar_sw',4) - output_sice_f(OutputFolder+"bba.dat",'albedo_bb_spherical_sw',5) - output_sice_f(OutputFolder+"size.dat",'D',4) - output_sice_f(OutputFolder+"size.dat",'area',5) - output_sice_f(OutputFolder+"size.dat",'al',6) - output_sice_f(OutputFolder+"size.dat",'r0',7) - output_sice_f(OutputFolder+"bba_alex_reduced.dat",'diagnostic_retrieval',3) + ds_f.dropna('x', 'all').dropna('y', 'all').plot(ax=ax[0], vmin=vmin, vmax=vmax) + ax[0].set_title(code_ver1) + ds_p.dropna('x', 'all').dropna('y', 'all').plot(ax=ax[1], vmin=vmin, vmax=vmax) + ax[1].set_title(code_ver2) - for i in range(21): - output_sice_f(OutputFolder+"spherical_albedo.dat",'albedo_spectral_spherical_'+str(i+1).zfill(2),4+i) - output_sice_f(OutputFolder+"planar_albedo.dat",'albedo_spectral_planar_'+str(i+1).zfill(2),4+i) - output_sice_f(OutputFolder+"boar.dat",'rBRR_'+str(i+1),4+i) + (ds_p-ds_f).dropna('x', 'all').dropna('y', 'all').plot(ax=ax[2]) + ax[2].set_title(code_ver2+' - '+code_ver1) + + ax[3].plot(ds_f.values.flatten(), + ds_p.values.flatten(), + marker ='.', linestyle='None') + ax[3].plot([ds_f.min(), ds_f.max()], [ds_f.min(), ds_f.max()],color='k') + ax[3].set_xlabel(code_ver1) + ax[3].set_ylabel(code_ver2) - #%% plotting oulooks - try: - os.mkdir(OutputFolder+'plots') - except: - print('folder exist') - - import matplotlib.pyplot as plt - - # fig,ax=bl.heatmap_discrete(rio.open(OutputFolder+'diagnostic_retrieval.tif').read(1), - # 'diagnostic_retrieval ') - # ax.set_title(OutputFolder) - # fig.savefig(OutputFolder+'plots/diagnostic_retrieval.png',bbox_inches='tight') + - var_list = ('albedo_bb_planar_sw','albedo_bb_spherical_sw') - for i in range(len(var_list)): - var_1 = rio.open(OutputFolder+var_list[i]+'.tif').read(1) - plt.figure(figsize=(10,15)) - bl.heatmap(var_1,var_list[i], col_lim=(0, 1) ,cmap_in='jet') - plt.title(OutputFolder) - plt.savefig(OutputFolder+'plots/'+var_list[i]+'_diff.png',bbox_inches='tight') - plt.ioff() - for i in np.append(np.arange(11), np.arange(21)): - var_name = 'albedo_spectral_spherical_'+ str(i+1).zfill(2) - var_1 = rio.open(OutputFolder+var_name+'.tif').read(1) - plt.figure(figsize=(10,15)) - bl.heatmap(var_1,var_name, col_lim=(0, 1) ,cmap_in='jet') - plt.title(OutputFolder) - plt.savefig(OutputFolder+'plots/'+var_name+'.png',bbox_inches='tight') - plt.ion() \ No newline at end of file diff --git a/sice_io.py b/sice_io.py index 3d7cbfa..5ecdc71 100644 --- a/sice_io.py +++ b/sice_io.py @@ -267,34 +267,27 @@ def to_zarr(self, append_dim=None): output_path, _ = os.path.splitext(output_path) ds.to_zarr(output_path + '.OUT.zarr', mode=mode, append_dim=append_dim, encoding=encodings, consolidated=True) - # def to_csv(self): - # # %% Output - # print('\nText file output') - # # data_in = pd.read_csv(sys.argv[1]) - # data_out = data_in - # data_out['grain_diameter'] = self.diameter - # # data_out['snow_specific_area']=area - # data_out['al'] = self.al - # data_out['r0'] = self.r0 - # data_out['diagnostic_retrieval'] = self.isnow - # data_out['conc'] = self.conc - # data_out['albedo_bb_planar_sw'] = self.rp3 - # data_out['albedo_bb_spherical_sw'] = self.rs3 - # for i in np.append(np.arange(11), np.arange(15,21)): - # # for i in np.arange(21): - # data_out['albedo_spectral_spherical_' + str(i + 1).zfill(2)] = self.alb_sph[i,:,:] - # for i in np.append(np.arange(11), np.arange(15,21)): - # data_out['rBRR_'+str(i+1).zfill(2)] = self.rp[i,:,:] - # basename, ext = os.path.splitext(self.filename) - # data_out.to_csv(basename + '_out.csv') - def write_output(snow, OutputFolder): - file_name_list = {'BXXX': 'O3_SICE', 'diameter': 'grain_diameter', 'area': 'snow_specific_area', - 'al': 'al', 'r0': 'r0', 'isnow': 'diagnostic_retrieval', 'conc': 'conc', - 'rp3': 'albedo_bb_planar_sw', 'rs3': 'albedo_bb_spherical_sw'} - for var in ['diameter', 'area', 'rp3', 'rs3', 'isnow', 'r0', 'al']: - snow[var].unstack(dim='xy').transpose('y', 'x').rio.to_raster(os.path.join(OutputFolder, file_name_list[var] + '.tif')) - + file_name_list = { + "BXXX": "O3_SICE", + "diameter": "grain_diameter", + "area": "snow_specific_area", + "al": "al", + "r0": "r0", + "isnow": "isnow", + "conc": "conc", + "rp3": "albedo_bb_planar_sw", + "rs3": "albedo_bb_spherical_sw", + "factor": "factor", + } + print('Printing out:') + for var in ["diameter", "area", "rp3", "rs3", "isnow", "r0", "al"]: + print(var) + snow[var].unstack(dim="xy").transpose("y", "x").rio.to_raster( + os.path.join(OutputFolder, file_name_list[var] + ".tif") + ) + snow.alb_sph.sel(band=0).unstack(dim="xy").transpose("y", "x").rio.to_raster(OutputFolder+'/alb_sph_01_solved.tif') + snow.rp.sel(band=0).unstack(dim="xy").transpose("y", "x").rio.to_raster(OutputFolder+'/alb_pl_01_solved.tif') # for var in ['BXXX', ]: # var = OLCI_scene[var].unstack(dim='xy').transpose('y', 'x').rio.to_raster(os.path.join(OutputFolder, file_name_list[var] + '.tif')) diff --git a/sice_lib.py b/sice_lib.py index 180052d..2bbd587 100644 --- a/sice_lib.py +++ b/sice_lib.py @@ -118,37 +118,25 @@ def process(OLCI_scene, compute_polluted=True, **kwargs): def process_by_chunk(OLCI_scene, chunk_size=150000, compute_polluted=True): size = OLCI_scene.sza.shape[0] nchunks = int(max(np.floor(size / chunk_size), 1)) - OLCI_chunks = OLCI_scene.chunk({'band': 21, 'xy': chunk_size}) + OLCI_chunks = OLCI_scene.chunk({"band": 21, "xy": chunk_size}) # snow_chunks = OLCI_chunks.map_blocks(process,kwargs={}, template = snow_template) - xy_chunk_indexes = np.array(OLCI_chunks.chunks['xy']).cumsum() - - diameter = [] - area = [] - rp3 = [] - rs3 = [] - isnow = [] - r0 = [] - al = [] - for i in range(len(xy_chunk_indexes)-1): + xy_chunk_indexes = np.array(OLCI_chunks.chunks["xy"]).cumsum() + + snow = xr.Dataset() + + for i in range(len(xy_chunk_indexes) - 1): print(f"{i+1} / {nchunks}") - chunk = OLCI_scene.isel(xy=slice(xy_chunk_indexes[i], xy_chunk_indexes[i+1])) + # define chunk + chunk = OLCI_scene.isel(xy=slice(xy_chunk_indexes[i], xy_chunk_indexes[i + 1])) + # process chunk snow_chunk = process(chunk) - diameter.append(snow_chunk.diameter) - area.append(snow_chunk.area) - rp3.append(snow_chunk.rp3) - rs3.append(snow_chunk.rs3) - isnow.append(snow_chunk.isnow) - r0.append(snow_chunk.r0) - al.append(snow_chunk.r0) + + if i == 0: + snow = snow_chunk.copy() + else: + snow = xr.concat([snow, snow_chunk], dim="xy") del snow_chunk - snow = xr.Dataset() - snow['diameter'] = xr.concat(diameter, dim='xy') - snow['area'] = xr.concat(area, dim='xy') - snow['rp3'] = xr.concat(rp3, dim='xy') - snow['rs3'] = xr.concat(rs3, dim='xy') - snow['isnow'] = xr.concat(isnow, dim='xy') - snow['r0'] = xr.concat(r0, dim='xy') - snow['al'] = xr.concat(al, dim='xy') + return snow