# coding: utf8 #***************************************************** # PROGRAMME : GPP autocorrélation # AUTEUR : C. MAGAND adapted from B. DELORME # CREATION : 18/07/2016 # PYTHON : PYTHON2.7 # OBJECTIF : Comparaison des GPP - simu - Mjung #****************************************************** #**** Packages loading import netCDF4 import scipy import numpy as np import numpy.ma as ma import datetime as dt import calendar import pandas as pd import pratic as p from scipy.interpolate import griddata # zones data zones=['Boreal_North_America','Temperate_North_America','Europe','Boreal_Asia','Temperate_Asia','Tropical_Asia','Australia_and_NewZealand','North_Africa','South_Africa','Tropical_South_America','South_America_Temperate','global'] #zones=['Boreal North America','Temperate North America','Europe','Boreal Asia','Temperate Asia','Tropical Asia','Australia + NewZealand','North Africa','South Africa','Tropical South America','South America Temperate','global'] #ind_zones=[1,2,3,4,5,6,7,8,9,10,11,16] ind_zones=[0,1,2,3,4,5,6,7,8,9,10,15] nb_zones=len(ind_zones) #12 #*** Input files and directories #indir='/Users/cmagand/pro/PC/autocor/gpp/' indir='/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/' #obs_file=indir+'MJUNGOctober2014_gpp_flux_1990_2009_96x73.nc' obs_file=indir+'fco2_JENA-s81-V3-6_Oct2014-ext3_1981-2000_monthlymean_XYT_144x143_ymonmean_units.nc' #mod_file=indir+'NPv5.70PDctrl01_19530101_19721231_1M_nbp_time_TS.nc' mod_file1=indir+'IPSL-CM5A-MR_19860101_20051231_1M_nbp_time_TS_units.nc' mod_file2=indir+'CM605-LR-pdCtrl-01_20600101_20791231_1M_nbp_time_TS_units.nc' mod_file3=indir+'CM605.calv-LR-pdCtrl-02_20900101_21091231_1M_nbp_time_TS_units.nc' mod_file4=indir+'CM605.dt20-LR-pdCtrl-02_20500101_20691231_1M_nbp_time_TS_units.nc' mod_file5=indir+'CM605.GUST-LR-pdCtrl-01_20800101_20991231_1M_nbp_time_TS_units.nc' mod_file6=indir+'CM605.NOSU-LR-pdCtrl-03_19800101_19991231_1M_nbp_time_TS_units.nc' mod_file7=indir+'CM605.THC1-LR-pdCtrl-01_21000101_21191231_1M_nbp_time_TS_units.nc' mod_file8=indir+'CM605.Z0-LR-pdCtrl-01_20800101_20991231_1M_nbp_time_TS_units.nc' mod_file9=indir+'NPv5.70PDctrl01_19530101_19721231_1M_nbp_time_TS_units.nc' mod_file10=indir+'v5.71vd2qsp375_21020101_21211231_1M_nbp_time_TS_units.nc' area_file=indir+'CM605.calv-LR-pdCtrl-02_20900101_21091231_1M_nbp_time.nc' #mod_file11=indir+'FG1low.Choi.CN3238.RUN.NewX_19940101_20131231_1M_nbp_time_TS_units.nc' #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/IPSL-CM5A-MR_19860101_20051231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/CM605-LR-pdCtrl-01_20600101_20791231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/CM605.calv-LR-pdCtrl-02_20900101_21091231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/CM605.dt20-LR-pdCtrl-02_20500101_20691231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/CM605.GUST-LR-pdCtrl-01_20800101_20991231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/CM605.NOSU-LR-pdCtrl-03_19800101_19991231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/CM605.THC1-LR-pdCtrl-01_21000101_21191231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/CM605.Z0-LR-pdCtrl-01_20800101_20991231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/NPv5.70PDctrl01_19530101_19721231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/v5.71vd2qsp375_21020101_21211231_1M_nbp_time_TS.nc #/data/cpipsl/CMIP6/INPUT/nbp/CMIP6-IPSL/SEMI-VARIOGRAM/FG1low.Choi.CN3238.RUN.NewX_19940101_20131231_1M_nbp_time_TS.nc #mod_file=indir+'v5.67PDay01_19820101_20011231_1M_gpp_time_96x73.nc' #mask_05 = "/Users/cmagand/pro/CMIP5/TOOLS/regions_mask_extend_land_96x73.nc" mask_05=indir+'regions_mask_extend_land_oct2015_144x143.nc' # get data # Coordinates nc = netCDF4.Dataset(obs_file,"r") lons = nc.variables["lon"][:]-180 lats = nc.variables["lat"][:] lon_ = len(nc.dimensions["lon"]) lat_ = len(nc.dimensions["lat"]) map_obs=nc.variables["nbp_units"][:] #*0.001/12. #*86400. # nc.close() #print map_obs[10][10] # Areas of cells nc= netCDF4.Dataset(area_file,"r") areas=nc.variables["Areas"][:] frac=nc.variables["CONTFRAC"][:] areacell=areas*frac areacell=areacell.filled(0) #map_mod1=nc.variables["nbp_units"][:]*86400.*(365./12.)*1000. #*86400.*365./1000. nc.close() nc1= netCDF4.Dataset(mod_file1,"r") #areas=nc1.variables["Areas"][:] #frac=nc1.variables["CONTFRAC"][:] map_mod1=nc1.variables["nbp_units"][:]*(-1.) #*areacell #*86400.*(365./12.)*1000. #*86400.*365./1000. nc2= netCDF4.Dataset(mod_file2,"r") #areas=nc2.variables["Areas"][:] #frac=nc2.variables["CONTFRAC"][:] map_mod2=nc2.variables["nbp_units"][:] #*areacell #*86400.*(365./12.)*1000. #*86400.*365./1000. nc3= netCDF4.Dataset(mod_file3,"r") map_mod3=nc3.variables["nbp_units"][:] #*areacell #*86400.*(365./12.)*1000. #*86400.*365./1000. nc4= netCDF4.Dataset(mod_file4,"r") map_mod4=nc4.variables["nbp_units"][:] #*areacell #*86400.*(365./12.)*1000. #*86400.*365./1000. nc5= netCDF4.Dataset(mod_file5,"r") map_mod5=nc5.variables["nbp_units"][:] #*areacell #*86400.*(365./12.)*1000. #*86400.*365./1000. nc6= netCDF4.Dataset(mod_file6,"r") map_mod6=nc6.variables["nbp_units"][:] #*areacell #*86400.*(365./12.)*1000. #*86400.*365./1000. nc7= netCDF4.Dataset(mod_file7,"r") map_mod7=nc7.variables["nbp_units"][:] #*areacell # *86400.*(365./12.)*1000. #*86400.*365./1000. nc8= netCDF4.Dataset(mod_file8,"r") map_mod8=nc8.variables["nbp_units"][:] #*areacell # *86400.*(365./12.)*1000. #*86400.*365./1000. nc9= netCDF4.Dataset(mod_file9,"r") map_mod9=nc9.variables["nbp_units"][:] #*areacell #*(-1.)*86400.*(365./12.)*1000. #*86400.*365./1000. nc10= netCDF4.Dataset(mod_file10,"r") map_mod10=nc10.variables["nbp_units"][:] #*areacell #*(-1.)*86400.*(365./12.)*1000. #*86400.*365./1000. #nc11= netCDF4.Dataset(mod_file11,"r") #map_mod11=nc11.variables["nbp_units"][:]*(-1.)*86400.*(365./12.)*1000. #*86400.*365./1000. #x = ma.array(data = map_mod, mask = map_mod_mask, fill_value = 0) # for i in range(lon_): # for j in range(lat_): # #print i # #print j # #print (gpp_obs_m[i,j]) # if x[i,j]: #map_mod.filled(9999) #print map_mod[10][10] #nc.close() nc1.close() nc2.close() nc3.close() nc4.close() nc5.close() nc6.close() nc7.close() nc8.close() nc9.close() nc10.close() #nc11.close() #areacell=areas*frac #areacell=areacell.filled(0) nmonths=map_obs.shape[0] print "additions" map_mod1[map_mod1>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod1=map_mod1+map_obs-map_obs map_obs=map_obs+map_mod1-map_mod1 #print map_obs[10][10] #map_mod1_mask = map_mod1.mask ##print map_mod_mask #map_mod1=map_mod1.filled(np.NaN) #map_obs_mask = map_obs.mask ##print map_mod_mask #map_obs=map_obs.filled(np.NaN) #print map_mod1[10][10] map_mod2[map_mod2>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod2=map_mod2+map_obs-map_obs map_obs=map_obs+map_mod2-map_mod2 #print map_obs[10][10] map_mod2_mask = map_mod2.mask #print map_mod_mask map_mod2=map_mod2.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod_mask map_obs=map_obs.filled(np.NaN) #print map_mod2[10][10] map_mod3[map_mod3>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod3=map_mod3+map_obs-map_obs map_obs=map_obs+map_mod3-map_mod3 #print map_obs[10][10] map_mod3_mask = map_mod3.mask #print map_mod_mask map_mod3=map_mod3.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod_mask map_obs=map_obs.filled(np.NaN) #print map_mod3[10][10] map_mod4[map_mod4>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod4=map_mod4+map_obs-map_obs map_obs=map_obs+map_mod4-map_mod4 #print map_obs[10][10] map_mod4_mask = map_mod4.mask #print map_mod_mask map_mod4=map_mod4.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod4_mask map_obs=map_obs.filled(np.NaN) #print map_mod4[10][10] map_mod5[map_mod5>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod5=map_mod5+map_obs-map_obs map_obs=map_obs+map_mod5-map_mod5 #print map_obs[10][10] map_mod5_mask = map_mod5.mask #print map_mod_mask map_mod=map_mod5.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod5_mask map_obs=map_obs.filled(np.NaN) #print map_mod5[10][10] map_mod6[map_mod6>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod6=map_mod6+map_obs-map_obs map_obs=map_obs+map_mod6-map_mod6 #print map_obs[10][10] map_mod6_mask = map_mod6.mask #print map_mod_mask map_mod6=map_mod6.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod6_mask map_obs=map_obs.filled(np.NaN) #print map_mod6[10][10] map_mod7[map_mod7>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod7=map_mod7+map_obs-map_obs map_obs=map_obs+map_mod7-map_mod7 #print map_obs[10][10] map_mod7_mask = map_mod7.mask #print map_mod_mask map_mod7=map_mod7.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod7_mask map_obs=map_obs.filled(np.NaN) #print map_mod7[10][10] map_mod8[map_mod8>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod8=map_mod8+map_obs-map_obs map_obs=map_obs+map_mod8-map_mod8 #print map_obs[10][10] map_mod8_mask = map_mod8.mask #print map_mod_mask map_mod8=map_mod8.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod8_mask map_obs=map_obs.filled(np.NaN) #print map_mod8[10][10] map_mod9[map_mod9>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod9=map_mod9+map_obs-map_obs map_obs=map_obs+map_mod9-map_mod9 #print map_obs[10][10] map_mod9_mask = map_mod9.mask #print map_mod_mask map_mod9=map_mod9.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod9_mask map_obs=map_obs.filled(np.NaN) #print map_mod9[10][10] map_mod10[map_mod10>1e30]=np.nan map_obs[map_obs<-1e30]=np.nan map_mod10=map_mod10+map_obs-map_obs map_obs=map_obs+map_mod10-map_mod10 #print map_obs[10][10] map_mod10_mask = map_mod10.mask #print map_mod_mask map_mod10=map_mod10.filled(np.NaN) map_obs_mask = map_obs.mask #print map_mod10_mask map_obs=map_obs.filled(np.NaN) #print map_mod10[10][10] #map_mod11[map_mod11>1e30]=np.nan #map_obs[map_obs<-1e30]=np.nan #map_mod11=map_mod11+map_obs-map_obs #map_obs=map_obs+map_mod11-map_mod11 ##print map_obs[10][10] #map_mod11_mask = map_mod11.mask ##print map_mod_mask #map_mod11=map_mod11.filled(np.NaN) #map_obs_mask = map_obs.mask ##print map_mod11_mask #map_obs=map_obs.filled(np.NaN) #print map_mod11[10][10] # Temporal Autocorrelation cut=7 print "temporal autocorrelation" #p.autocor(map_obs,map_mod,mask_05,nb_zones,zones,'FD',ind_zones,cut,nmonths) # Spatial autocorrelation days=[15,45,75,105,135,165,195,225,255,285,315,345] nlag=30 lmax=3000 print "semivario" #print map_obs #p.semivario(map_obs,map_mod,mask_05,nb_zones,zones,'NBP',ind_zones,days,nlag,lmax) p.semivario_multi(map_obs,map_mod1,map_mod2,map_mod3,map_mod4,map_mod5,map_mod6,map_mod7,map_mod8,map_mod9,map_mod10,mask_05,nb_zones,zones,'NBP',ind_zones,days,nlag,lmax)