From 7db5effa171502775245e0cc9ebf03a063200d9b Mon Sep 17 00:00:00 2001 From: sherimickelson Date: Mon, 2 Jul 2018 20:15:10 -0600 Subject: [PATCH 1/4] Add more user functions Add code to iconform to prevent it from creating input specs for dim variables --- scripts/iconform | 16 +- source/pyconform/modules/commonfunctions.py | 253 ++++++++++++++++++++ 2 files changed, 262 insertions(+), 7 deletions(-) diff --git a/scripts/iconform b/scripts/iconform index 5fdc7cfe..46a0d651 100755 --- a/scripts/iconform +++ b/scripts/iconform @@ -514,15 +514,17 @@ def create_output(exp_dict, definitions, input_glob, attributes, output_path, ar with open(f, 'w') as outfile: json.dump(t, outfile, sort_keys=True, indent=4) else: + ignore = ['latitude','longitude','olevel','plev19','time','time1','alevhalf','ygre','xgre','vegtype','spectband','areacellg','alevel','xant','yant','rho','tau','plev3','gridlatitude','plev39','plev4','plev27','plev3','plev7h','plev7c','plev8','plev7h'] for n,t in TableSpec.iteritems(): for vn,var in t.iteritems(): - varD={} - varD[vn]=var - for d in var["dimensions"]: - varD[d]=t[d] - f = output_path + "/" + experiment + '_' + n + '_' + vn + '_spec.json' - with open(f, 'w') as outfile: - json.dump(varD, outfile, sort_keys=True, indent=4) + if vn not in ignore: + varD={} + varD[vn]=var + for d in var["dimensions"]: + varD[d]=t[d] + f = output_path + "/" + experiment + '_' + n + '_' + vn + '_spec.json' + with open(f, 'w') as outfile: + json.dump(varD, outfile, sort_keys=True, indent=4) #f1 = output_path + '/MISSING_DEFS.json' #with open(f1, 'w') as outfile: diff --git a/source/pyconform/modules/commonfunctions.py b/source/pyconform/modules/commonfunctions.py index 029d3566..d89898ef 100644 --- a/source/pyconform/modules/commonfunctions.py +++ b/source/pyconform/modules/commonfunctions.py @@ -259,6 +259,74 @@ def __getitem__(self, index): return PhysArray(a, name = new_name, units=p_data.units) + +#=================================================================================================== +# sftofFunction +#=================================================================================================== +class sftofFunction(Function): + key = 'sftof' + + def __init__(self, KMT): + super(sftofFunction, self).__init__(KMT) + + def __getitem__(self, index): + p_KMT = self.arguments[0][index] + + if index is None: + return PhysArray(np.zeros((0,0)), dimensions=[p_KMT.dimensions[0],p_KMT.dimensions[1]]) + + KMT = p_KMT.data + + a = np.zeros((KMT.shape[0],KMT.shape[1])) + + for j in range(KMT.shape[0]): + for i in range(KMT.shape[1]): + if KMT[j,i] > 0: + a[j,i] = 1 + + new_name = 'sftof({})'.format( p_KMT.name) + + return PhysArray(a, name = new_name, dimensions=[p_KMT.dimensions[0],p_KMT.dimensions[1]], units=p_KMT.units) + + + +#=================================================================================================== +# POP_bottom_layer_multFunction +#=================================================================================================== +class POP_bottom_layer_multaddFunction(Function): + key = 'POP_bottom_layer_multadd' + + def __init__(self, KMT, data1, data2): + super(POP_bottom_layer_multaddFunction, self).__init__(KMT, data1, data2) + + def __getitem__(self, index): + p_KMT = self.arguments[0][index] + p_data1 = self.arguments[1][index] + p_data2 = self.arguments[2][index] + + data1 = p_data1.data + data2 = p_data2.data + KMT = p_KMT.data + + a1 = np.zeros((p_data2.shape[0],p_data2.shape[2],p_data2.shape[3])) + a2 = np.zeros((p_data2.shape[0])) + + for t in range(p_data2.shape[0]): + for j in range(KMT.shape[0]): + for i in range(KMT.shape[1]): + k = KMT[j,i]-1 + if data2[t,k,j,i] < 1e+16: + a1[t,j,i] = data1[k]*data2[t,k,j,i] + #print a1[t,j,i],data1[k],data2[t,k,j,i] + for t in range(p_data2.shape[0]): + a2[t] = np.ma.sum(a1[t,:,:]) + #print a2[t] + + return PhysArray(a2, dimensions=[p_data2.dimensions[0]], units=p_data2.units) + + + + #=================================================================================================== # masked_invalidFunction #=================================================================================================== @@ -357,4 +425,189 @@ def __getitem__(self, index): return PhysArray(a, dimensions=[a1.dimensions[0],a1.dimensions[1],a1.dimensions[2]], units=var.units) +#=================================================================================================== +# cice_regions +#=================================================================================================== +class cice_regionsFunction(Function): + key = 'cice_regions' + + def __init__(self, p_aice, p_uvel, p_vvel, p_HTE, p_HTN, p_siline, multiple): + super(cice_regionsFunction, self).__init__(p_aice, p_uvel, p_vvel, p_HTE, p_HTN, p_siline, multiple) + + def __getitem__(self, index): + p_aice = self.arguments[0][index] + p_uvel = self.arguments[1][index] + p_vvel = self.arguments[2][index] + p_HTE = self.arguments[3][index] + p_HTN = self.arguments[4][index] + p_siline = self.arguments[5][index] + multiple = self.arguments[6] + + aice=p_aice + uvel=p_uvel + vvel=p_vvel + HTE=p_HTE + HTN=p_HTN + siline=p_siline + a = np.ma.zeros((aice.data.shape[0],siline.data.shape[0])) + + for t in range(aice.data.shape[0]): + #1 + i = 92 + for j in range(370, 381): + if aice[t,j,i] < 1e+16 and aice[t,j,i+1] < 1e+16 and HTE[j,i] < 1e+16 and uvel[t,j,i] < 1e+16 and uvel[t,j-1,i] < 1e+16: + a[t,0] += 0.5*(aice[t,j,i]+aice[t,j,i+1])*0.5*(HTE[j,i]*uvel[t,j,i]+HTE[j,i]*uvel[t,j-1,i]) + #2 + i = 214 + for j in range(375,377): + if aice[t,j,i] < 1e+16 and aice[t,j,i+1] < 1e+16 and HTE[j,i] < 1e+16 and uvel[t,j,i] < 1e+16 and uvel[t,j-1,i] < 1e+16 and aice[t,j+1,i] < 1e+16 and HTN[j,i] < 1e+16 and vvel[t,j,i] < 1e+16 and vvel[t,j,i-1] < 1e+16: + a[t,1] += 0.5*(aice[t,j,i]+aice[t,j,i+1])*0.5*(HTE[j,i]*uvel[t,j,i]+HTE[j,i]*uvel[t,j-1,i]) + 0.5*(aice[t,j,i]+aice[t,j+1,i])*0.5*(HTN[j,i]*vvel[t,j,i]+HTN[j,i]*vvel[t,j,i-1]) + j = 366 + for i in range(240,244): + if aice[t,j,i] < 1e+16 and aice[t,j,i+1] < 1e+16 and HTE[j,i] < 1e+16 and uvel[t,j,i] < 1e+16 and uvel[t,j-1,i] < 1e+16 and aice[t,j+1,i] < 1e+16 and HTN[j,i] < 1e+16 and vvel[t,j,i] < 1e+16 and vvel[t,j,i-1] < 1e+16: + a[t,1] += 0.5*(aice[t,j,i]+aice[t,j,i+1])*0.5*(HTE[j,i]*uvel[t,j,i]+HTE[j,i]*uvel[t,j-1,i]) + 0.5*(aice[t,j,i]+aice[t,j+1,i])*0.5*(HTN[j,i]*vvel[t,j,i]+HTN[j,i]*vvel[t,j,i-1]) + #3 + i = 85 + for j in range(344,366): + if aice[t,j,i] < 1e+16 and aice[t,j,i+1] < 1e+16 and HTE[j,i] < 1e+16 and uvel[t,j,i] < 1e+16 and uvel[t,j-1,i] < 1e+16: + a[t,2] += 0.5*(aice[t,j,i]+aice[t,j,i+1])*0.5*(HTE[j,i]*uvel[t,j,i]+HTE[j,i]*uvel[t,j-1,i]) + #4 + j = 333 + for j in range(198,201): + if aice[t,j,i] < 1e+16 and aice[t,j+1,i] < 1e+16 and HTN[j,i] < 1e+16 and vvel[t,j,i] < 1e+16 and vvel[t,j,i-1] < 1e+16: + a[t,3] += 0.5*(aice[t,j,i]+aice[t,j+1,i])*0.5*(HTN[j,i]*vvel[t,j,i]+HTN[j,i]*vvel[t,j,i-1]) + + a = a*multiple + + return PhysArray(a, dimensions=[p_aice.dimensions[0],p_siline.dimensions[0]], units=uvel.units) + + +#=================================================================================================== +# reduce_luFunction +#=================================================================================================== +class reduce_luFunction(Function): + key = 'reduce_lu' + + # np.where(x < 5, x, -1) + + def __init__(self, p_data,p_lu): + super(reduce_luFunction, self).__init__(p_data,p_lu) + + def __getitem__(self, index): + p_data = self.arguments[0][index] + p_lu = self.arguments[1][index] + + #if index is None: + # return PhysArray(p_data.data, dimensions=[p_data.dimensions[0],p_lu.dimensions[0],p_data.dimensions[2],p_data.dimensions[3]]) + + data = p_data.data + lu = p_lu.data + + data2 = np.ma.zeros((data.shape[0],4,data.shape[2],data.shape[3])) + + for t in range(data.shape[0]): + for x in range(data.shape[2]): + for y in range(data.shape[3]): + data2[t,0,x,y] = data[t,0,x,y] + data2[t,1,x,y] = 0 + data2[t,2,x,y] = data[t,1,x,y] + data2[t,3,x,y] = data[t,6,x,y]+data[t,7,x,y]+data[t,8,x,y] + data2[data2>=1e+16] = 1e+20 + + return PhysArray(data2, dimensions=[p_data.dimensions[0],p_lu.dimensions[0],p_data.dimensions[2],p_data.dimensions[3]], units=p_data.units) + + +#=================================================================================================== +# soilpoolsFunction +#=================================================================================================== +class soilpoolsFunction(Function): + key = 'soilpools' + + def __init__(self, p_data1,p_data2,p_data3,p_soilpool): + super(soilpoolsFunction, self).__init__(p_data1,p_data2,p_data3,p_soilpool) + + def __getitem__(self, index): + p_data1 = self.arguments[0][index] + p_data2 = self.arguments[1][index] + p_data3 = self.arguments[2][index] + p_soilpool = self.arguments[3][index] + + data1 = p_data1.data + data2 = p_data2.data + data3 = p_data3.data + soilpool = p_soilpool.data + + data = np.ma.zeros((data1.shape[0],3,data1.shape[1],data1.shape[2])) + + data[:,0,:,:] = data1 + data[:,1,:,:] = data2 + data[:,2,:,:] = data3 + + data[data>=1e+16] = 1e+20 + data = np.ma.masked_values(data, 1e+20) + + return PhysArray(data, dimensions=[p_data1.dimensions[0],p_soilpool.dimensions[0],p_data1.dimensions[1],p_data1.dimensions[2]], units=p_data1.units) + + +#=================================================================================================== +# expand_latlonFunction +#=================================================================================================== +class expand_latlonFunction(Function): + key = 'expand_latlon' + + def __init__(self, p_data1,p_lat,p_lon): + super(expand_latlonFunction, self).__init__(p_data1,p_lat,p_lon) + + def __getitem__(self, index): + p_data1 = self.arguments[0][index] + p_lat = self.arguments[1][index] + p_lon = self.arguments[2][index] + + data1 = p_data1.data + lat = p_lat.data + lon = p_lon.data + + data = np.ma.zeros((data1.shape[0],lat.shape[0],lon.shape[0])) + + for x in range(lat.shape[0]): + for y in range(lon.shape[0]): + data[:,x,y] = data1 + + data[data>=1e+16] = 1e+20 + data = np.ma.masked_values(data, 1e+20) + + return PhysArray(data, dimensions=[p_data1.dimensions[0],p_lat.dimensions[0],p_lon.dimensions[0]], units=p_data1.units) + + +#=================================================================================================== +# ocean_basinFunction +#=================================================================================================== +class ocean_basinFunction(Function): + key = 'ocean_basin' + + def __init__(self, p_data1, p_comp, p_basin): + super(ocean_basinFunction, self).__init__(p_data1, p_comp, p_basin) + + def __getitem__(self, index): + p_data1 = self.arguments[0][index] + p_comp = self.arguments[1] + p_basin = self.arguments[2][index] + + data1 = p_data1.data + comp = int(p_comp) + basin = p_basin.data + + data = np.ma.zeros((data1.shape[0],data1.shape[4],data1.shape[3],basin.shape[0])) + + for t in range(data1.shape[0]): + for x in range(data1.shape[4]): + for y in range(data1.shape[3]): + data[t,x,y,0] = data1[t,1,comp,y,x] + data[t,x,y,1] = data1[t,0,comp,y,x]-data1[t,1,comp,y,x] + data[t,x,y,2] = data1[t,0,comp,y,x] + + data[data>=1e+16] = 1e+20 + data = np.ma.masked_values(data, 1e+20) + + return PhysArray(data, dimensions=[p_data1.dimensions[0],p_data1.dimensions[4],p_data1.dimensions[3],p_basin.dimensions[0]], units=p_data1.units) From 18bdddad17feba690d32c7579aba620211bd6f93 Mon Sep 17 00:00:00 2001 From: sherimickelson Date: Mon, 9 Jul 2018 14:41:20 -0600 Subject: [PATCH 2/4] Add names to several functions in commonfunctions.py Catch for None entries within an input json file within the iconform script --- scripts/iconform | 4 ++++ source/pyconform/modules/commonfunctions.py | 25 ++++++++++++++------- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/scripts/iconform b/scripts/iconform index 46a0d651..a712677a 100755 --- a/scripts/iconform +++ b/scripts/iconform @@ -467,6 +467,7 @@ def create_output(exp_dict, definitions, input_glob, attributes, output_path, ar ts_key = None mip = d['mipTable'] if mip in definitions.keys(): + ig = "" if v in definitions[mip].keys(): if "N/A" in definitions[mip][v].upper(): v_def = "" @@ -614,6 +615,9 @@ def main(argv=None): if "json" in gaFile: with open(gaFile) as gaF: ga = json.load(gaF) + for k in ga.keys(): + if ga[k] is None: + ga[k] = "" attributes.update(ga) else: with open(gaFile) as y_attributes: diff --git a/source/pyconform/modules/commonfunctions.py b/source/pyconform/modules/commonfunctions.py index d89898ef..b9501b05 100644 --- a/source/pyconform/modules/commonfunctions.py +++ b/source/pyconform/modules/commonfunctions.py @@ -322,7 +322,8 @@ def __getitem__(self, index): a2[t] = np.ma.sum(a1[t,:,:]) #print a2[t] - return PhysArray(a2, dimensions=[p_data2.dimensions[0]], units=p_data2.units) + new_name = 'POP_bottom_layer_multadd({}{}{})'.format(KMT.name, data1.name, data2.name) + return PhysArray(a2, name=new_name, dimensions=[p_data2.dimensions[0]], units=p_data2.units) @@ -422,7 +423,9 @@ def __getitem__(self, index): a[t,:,:] = np.ma.where(a1[t,:,:] < a2, var, value) elif '>' in condition: a[t,:,:] = np.ma.where(a1[t,:,:] > a2, var, value) - return PhysArray(a, dimensions=[a1.dimensions[0],a1.dimensions[1],a1.dimensions[2]], units=var.units) + + new_name = 'cice_where()'.format() + return PhysArray(a, name=new_name, dimensions=[a1.dimensions[0],a1.dimensions[1],a1.dimensions[2]], units=var.units) #=================================================================================================== @@ -473,13 +476,14 @@ def __getitem__(self, index): a[t,2] += 0.5*(aice[t,j,i]+aice[t,j,i+1])*0.5*(HTE[j,i]*uvel[t,j,i]+HTE[j,i]*uvel[t,j-1,i]) #4 j = 333 - for j in range(198,201): + for i in range(198,201): if aice[t,j,i] < 1e+16 and aice[t,j+1,i] < 1e+16 and HTN[j,i] < 1e+16 and vvel[t,j,i] < 1e+16 and vvel[t,j,i-1] < 1e+16: a[t,3] += 0.5*(aice[t,j,i]+aice[t,j+1,i])*0.5*(HTN[j,i]*vvel[t,j,i]+HTN[j,i]*vvel[t,j,i-1]) a = a*multiple - return PhysArray(a, dimensions=[p_aice.dimensions[0],p_siline.dimensions[0]], units=uvel.units) + new_name = 'cice_regions()'.format() + return PhysArray(a, name=new_name, dimensions=[p_aice.dimensions[0],p_siline.dimensions[0]], units=uvel.units) #=================================================================================================== @@ -514,7 +518,8 @@ def __getitem__(self, index): data2[t,3,x,y] = data[t,6,x,y]+data[t,7,x,y]+data[t,8,x,y] data2[data2>=1e+16] = 1e+20 - return PhysArray(data2, dimensions=[p_data.dimensions[0],p_lu.dimensions[0],p_data.dimensions[2],p_data.dimensions[3]], units=p_data.units) + new_name = 'reduce_lu({}{})'.format(p_data.name,p_lu.name) + return PhysArray(data2, name=new_name, dimensions=[p_data.dimensions[0],p_lu.dimensions[0],p_data.dimensions[2],p_data.dimensions[3]], units=p_data.units) #=================================================================================================== @@ -546,7 +551,8 @@ def __getitem__(self, index): data[data>=1e+16] = 1e+20 data = np.ma.masked_values(data, 1e+20) - return PhysArray(data, dimensions=[p_data1.dimensions[0],p_soilpool.dimensions[0],p_data1.dimensions[1],p_data1.dimensions[2]], units=p_data1.units) + new_name = 'soilpools({}{}{}{})'.format(p_data1.name,p_data2.name,p_data3.name,p_soilpool.name) + return PhysArray(data, name=new_name, dimensions=[p_data1.dimensions[0],p_soilpool.dimensions[0],p_data1.dimensions[1],p_data1.dimensions[2]], units=p_data1.units) #=================================================================================================== @@ -576,7 +582,8 @@ def __getitem__(self, index): data[data>=1e+16] = 1e+20 data = np.ma.masked_values(data, 1e+20) - return PhysArray(data, dimensions=[p_data1.dimensions[0],p_lat.dimensions[0],p_lon.dimensions[0]], units=p_data1.units) + new_name = 'expand_latlon({}{}{})'.format( p_data1.name, p_lat.name, p_lon.name) + return PhysArray(data, name=new_name, dimensions=[p_data1.dimensions[0],p_lat.dimensions[0],p_lon.dimensions[0]], units=p_data1.units) #=================================================================================================== @@ -609,5 +616,7 @@ def __getitem__(self, index): data[data>=1e+16] = 1e+20 data = np.ma.masked_values(data, 1e+20) - return PhysArray(data, dimensions=[p_data1.dimensions[0],p_data1.dimensions[4],p_data1.dimensions[3],p_basin.dimensions[0]], units=p_data1.units) + new_name = 'ocean_basin({}{})'.format(p_data1.name, p_basin.name) + return PhysArray(data, name=new_name, dimensions=[p_data1.dimensions[0],p_data1.dimensions[4],p_data1.dimensions[3],p_basin.dimensions[0]], units=p_data1.units) + From 745334e8b28736f523b434c995ea0d0b33bd7d7e Mon Sep 17 00:00:00 2001 From: Kevin Paul Date: Thu, 12 Jul 2018 14:33:31 -0600 Subject: [PATCH 3/4] Semi-optimized CLM PFT to vegtype function (73% speedup) --- .../modules/CLM_pft_to_CMIP6_vegtype.py | 340 +++++++----------- 1 file changed, 136 insertions(+), 204 deletions(-) diff --git a/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py b/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py index e62e8f9d..759cc1ef 100644 --- a/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py +++ b/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py @@ -1,9 +1,8 @@ #! /usr/bin/env python -import time, sys import numpy as np -from pyconform.physarray import PhysArray, UnitsError, DimensionsError -from pyconform.functions import Function, is_constant +from pyconform.physarray import PhysArray +from pyconform.functions import Function class CLM_pft_to_CMIP6_vegtype_Function(Function): @@ -11,18 +10,17 @@ class CLM_pft_to_CMIP6_vegtype_Function(Function): numargs = 18 def __init__(self, GPP, vegType, time, lat, lon, grid1d_ixy, grid1d_jxy, grid1d_lon, - grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, - pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, - pfts1d_wtgcell, pfts1d_wtlunit): + grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, + pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, + pfts1d_wtgcell, pfts1d_wtlunit): super(CLM_pft_to_CMIP6_vegtype_Function, self).__init__(GPP, vegType, time, lat, lon, grid1d_ixy, grid1d_jxy, grid1d_lon, - grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, - pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, - pfts1d_wtgcell, pfts1d_wtlunit) + grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, + pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, + pfts1d_wtgcell, pfts1d_wtlunit) def __getitem__(self, index): - pGPP = self.arguments[0][index] # vegType = grass, shrub, or tree vegType = self.arguments[1] @@ -32,22 +30,22 @@ def __getitem__(self, index): pgrid1d_ixy = self.arguments[5][index] pgrid1d_jxy = self.arguments[6][index] pgrid1d_lon = self.arguments[7][index] - pgrid1d_lat = self.arguments[8][index] + pgrid1d_lat = self.arguments[8][index] pland1d_lon = self.arguments[9][index] pland1d_lat = self.arguments[10][index] pland1d_ityplunit = self.arguments[11][index] - ppfts1d_lon = self.arguments[12][index] + ppfts1d_lon = self.arguments[12][index] ppfts1d_lat = self.arguments[13][index] ppfts1d_active = self.arguments[14][index] ppfts1d_itype_veg = self.arguments[15][index] - ppfts1d_wtgcell = self.arguments[16][index] + ppfts1d_wtgcell = self.arguments[16][index] ppfts1d_wtlunit = self.arguments[17][index] if index is None: - return PhysArray(np.zeros((0,0,0)), dimensions=[ptime.dimensions[0],plat.dimensions[0],plon.dimensions[0]]) - + return PhysArray(np.zeros((0, 0, 0)), dimensions=[ptime.dimensions[0], plat.dimensions[0], plon.dimensions[0]]) + GPP = pGPP.data - time = ptime.data + time = ptime.data lat = plat.data lon = plon.data grid1d_ixy = pgrid1d_ixy.data @@ -63,194 +61,128 @@ def __getitem__(self, index): pfts1d_itype_veg = ppfts1d_itype_veg.data pfts1d_wtgcell = ppfts1d_wtgcell.data pfts1d_wtlunit = ppfts1d_wtlunit.data - - - # vegType = grass, shrub, or tree - - # Tolerance check for weights summing to 1 - eps = 1.e-5 - - # If 1, pft is active - active_pft = 1 - - # If 1, landunit is veg - veg_lunit = 1 - - # C3 arctic grass, - # C3 non-arctic grass, - # C4 grass - beg_grass_pfts = 12 - end_grass_pfts = 14 - - # broadleaf evergreen shrub - temperate, - # broadleaf deciduous shrub - temperate, - # broadleaf deciduous shrub - boreal - beg_shrub_pfts = 9 - end_shrub_pfts = 11 - - # needleleaf evergreen tree - temperate, - # needleleaf evergreen tree - boreal, - # needleleaf deciduous tree - boreal, - # broadleaf evergreen tree - tropical, - # broadleaf evergreen tree - temperate, - # broadleaf deciduous tree - tropical, - # broadleaf deciduous tree - temperate, - # broadleaf deciduous tree - boreal - beg_tree_pfts = 1 - end_tree_pfts = 8 - - # Will contain weighted average for grass pfts on 2d grid - varo_vegType = np.zeros([len(time),len(lat),len(lon)]) - tu = np.stack((pfts1d_lon,pfts1d_lat, pfts1d_active), axis=1) - - ind = np.stack((grid1d_ixy,grid1d_jxy), axis=1) - - lu = np.stack((land1d_lon,land1d_lat,land1d_ityplunit), axis=1) - - # Loop over lat/lons - for ixy in range(len(lon)): - for jxy in range(len(lat)): - grid_indx = -99 - # 1d grid index - ind_comp = (ixy+1,jxy+1) - gi = np.where(np.all(ind==ind_comp, axis=1))[0] - if len(gi) > 0: - grid_indx = gi[0] - - # Check for valid land gridcell - if grid_indx != -99: - - # Gridcell lat/lons - grid1d_lon_pt = grid1d_lon[grid_indx] - grid1d_lat_pt = grid1d_lat[grid_indx] - - # veg landunit index for this gridcell - t_var = (grid1d_lon_pt, grid1d_lat_pt, veg_lunit) - landunit_indx = np.where(np.all(t_var == lu, axis=1))[0] - - # Check for valid veg landunit - if landunit_indx.size > 0: - if 'grass' in vegType: - t_var = (grid1d_lon_pt,grid1d_lat_pt, active_pft) - pft_indx = np.where( np.all(t_var == tu, axis=1) * (pfts1d_wtgcell > 0.) * (pfts1d_itype_veg >= beg_grass_pfts) * (pfts1d_itype_veg <= end_grass_pfts))[0] - elif 'shrub' in vegType: - t_var = (grid1d_lon_pt,grid1d_lat_pt, active_pft) - pft_indx = np.where( np.all(t_var==tu, axis=1) * (pfts1d_wtgcell > 0.) * (pfts1d_itype_veg >= beg_shrub_pfts) * (pfts1d_itype_veg <= end_shrub_pfts))[0] - elif 'tree' in vegType: - t_var = (grid1d_lon_pt,grid1d_lat_pt, active_pft) - pft_indx = np.where( np.all(t_var==tu, axis=1) * (pfts1d_wtgcell > 0.) * (pfts1d_itype_veg >= beg_tree_pfts) * (pfts1d_itype_veg <= end_tree_pfts))[0] - - # Check for valid pfts and compute weighted average - if pft_indx.size > 0: - if 'grass' in vegType: - pfts1d_wtlunit_grass = (pfts1d_wtlunit[pft_indx]).astype(np.float32) - dum = GPP[:,pft_indx] - weights = pfts1d_wtlunit_grass / np.sum(pfts1d_wtlunit_grass) - if np.absolute(1.-np.sum(weights)) > eps: - print("Weights do not sum to 1, exiting") - sys.exit(-1) - varo_vegType[:,jxy,ixy] = np.sum(dum * weights) - - elif 'shrub' in vegType: - pfts1d_wtlunit_shrub = (pfts1d_wtlunit[pft_indx]).astype(np.float32) - dum = GPP[:,pft_indx] - weights = pfts1d_wtlunit_shrub / np.sum(pfts1d_wtlunit_shrub) - varo_vegType[:,jxy,ixy] = np.sum(dum * weights) - - elif 'tree' in vegType: - pfts1d_wtlunit_tree = (pfts1d_wtlunit[pft_indx]).astype(np.float32) - dum = GPP[:,pft_indx] - weights = pfts1d_wtlunit_tree / np.sum(pfts1d_wtlunit_tree) - varo_vegType[:,jxy,ixy] = np.sum(dum * weights) - - else: - varo_vegType[:,jxy,ixy] = 1e+20 - else: - varo_vegType[:,jxy,ixy] = 1e+20 - else: - varo_vegType[:,jxy,ixy] = 1e+20 - - new_name = 'CLM_pft_to_CMIP6_vegtype({}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{})'.format(pGPP.name, - vegType, ptime.name, plat.name, plon.name, pgrid1d_ixy.name, pgrid1d_jxy.name, pgrid1d_lon.name, - pgrid1d_lat.name, pland1d_lon.name, pland1d_lat.name, pland1d_ityplunit.name, - ppfts1d_lon.name, ppfts1d_lat.name, ppfts1d_active.name, ppfts1d_itype_veg.name, - ppfts1d_wtgcell.name, ppfts1d_wtlunit.name) - - print 'FINISHED FUNCTION' - - varo_vegType[varo_vegType>=1e+16] = 1e+20 + # If 1, pft is active + active_pft = 1 + + # If 1, landunit is veg + veg_lunit = 1 + + # Veg-types and valid pfts + if 'tree' in vegType: + beg_pfts = 1 + end_pfts = 8 + elif 'shrub' in vegType: + beg_pfts = 9 + end_pfts = 11 + elif 'grass' in vegType: + beg_pfts = 12 + end_pfts = 14 + valid_pfts = ((pfts1d_wtgcell > 0.) * + (pfts1d_itype_veg >= beg_pfts) * + (pfts1d_itype_veg <= end_pfts)) + + # Will contain weighted average for grass pfts on 2d grid + varo_vegType = np.ones([len(time), len(lat), len(lon)]) * 1e20 + + vegtyp_lunit_indices = np.where(land1d_ityplunit == veg_lunit) + active_ptype_indices = np.where(pfts1d_active == active_pft) + + gcell_pts = np.stack((grid1d_lon, grid1d_lat), axis=1) + lunit_pts = np.stack( + (land1d_lon[vegtyp_lunit_indices], land1d_lat[vegtyp_lunit_indices]), axis=1) + pfts_pts = np.stack( + (pfts1d_lon[active_ptype_indices], pfts1d_lat[active_ptype_indices]), axis=1) + + gcell_indices = np.where(np.all(np.isin(gcell_pts, lunit_pts), axis=1)) + gcell_pts = gcell_pts[gcell_indices] + gcell_indices = np.where(np.all(np.isin(gcell_pts, pfts_pts), axis=1)) + + gcell_lon = grid1d_lon[gcell_indices] + gcell_lat = grid1d_lat[gcell_indices] + gcell_ixy = grid1d_ixy[gcell_indices] + gcell_jxy = grid1d_jxy[gcell_indices] + + for lon, lat, i, j in np.nditer([gcell_lon, gcell_lat, gcell_ixy, gcell_jxy]): + pft_idx = np.where( + np.all((lon, lat) == pfts_pts, axis=1) * valid_pfts)[0] + dum = GPP[:, pft_idx] + pfts1d_wtlunit_veg = (pfts1d_wtlunit[pft_idx]).astype(np.float32) + weights = pfts1d_wtlunit_veg / np.sum(pfts1d_wtlunit_veg) + varo_vegType[:, j - 1, i - 1] = np.sum(dum * weights) + + new_name = 'CLM_pft_to_CMIP6_vegtype({},{})'.format(pGPP.name, vegType) + + varo_vegType[varo_vegType >= 1e+16] = 1e+20 ma_varo_vegType = np.ma.masked_values(varo_vegType, 1e+20) - return PhysArray(ma_varo_vegType, name=new_name, units=pGPP.units) - - -def main(argv=None): - - import netCDF4 as nc - - sim = "clm50_r243_1deg_GSWP3V2_cropopt_nsc_emergeV2F_dailyo_hist" - f_in = sim+".clm2.h1.2005-01.nc" - f_out = sim+".clm2.h1veg.0001-01.nc" - f_dir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/DATA/" - f_outfir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/new/OUTDIR/" - - cdf_file = nc.Dataset(f_dir+f_in,"r") - - ntim = cdf_file.variables['time'][:] - nlat = cdf_file.variables['lat'][:] - nlon = cdf_file.variables['lon'][:] - - grid1d_ixy = cdf_file.variables['grid1d_ixy'][:] - grid1d_jxy = cdf_file.variables['grid1d_jxy'][:] - grid1d_lon = cdf_file.variables['grid1d_lon'][:] - grid1d_lat = cdf_file.variables['grid1d_lat'][:] - land1d_lon = cdf_file.variables['land1d_lon'][:] - land1d_lat = cdf_file.variables['land1d_lat'][:] - land1d_ityplunit = cdf_file.variables['land1d_ityplunit'][:] - pfts1d_lon = cdf_file.variables['pfts1d_lon'][:] - pfts1d_lat = cdf_file.variables['pfts1d_lat'][:] - pfts1d_active = cdf_file.variables['pfts1d_active'][:] - pfts1d_itype_veg = cdf_file.variables['pfts1d_itype_veg'][:] - pfts1d_wtgcell = cdf_file.variables['pfts1d_wtgcell'][:] - pfts1d_wtlunit = cdf_file.variables['pfts1d_wtlunit'][:] - - - GPP = cdf_file.variables['GPP'][:] - - cdf_file.close() - - out_file = nc.Dataset(f_outfir+f_out,"w") - - time = out_file.createDimension('time',None) - lat = out_file.createDimension('lat',len(nlat)) - lon = out_file.createDimension('lon',len(nlon)) - gppGrass = out_file.createVariable('gppGrass', 'f4', ('time', 'lat', 'lon'),fill_value=1.e36) - gppShrub = out_file.createVariable('gppShrub', 'f4', ('time', 'lat', 'lon'),fill_value=1.e36) - gppTree = out_file.createVariable('gppTree', 'f4', ('time', 'lat', 'lon'),fill_value=1.e36) - - print 'Looking for grass' - gppGrass[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'grass', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, - grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, - pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, - pfts1d_wtgcell, pfts1d_wtlunit) - print 'Looking for shrubs' - gppShrub[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'shrub', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, - grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, - pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, - pfts1d_wtgcell, pfts1d_wtlunit) - print 'Looking for trees' - gppTree[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'tree', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, - grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, - pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, - pfts1d_wtgcell, pfts1d_wtlunit) - - - out_file.close() - -if __name__ == '__main__': - main() - - - - + return PhysArray(ma_varo_vegType, name=new_name, units=pGPP.units) + + +# def main(): +# +# import netCDF4 as nc +# +# sim = "clm50_r243_1deg_GSWP3V2_cropopt_nsc_emergeV2F_dailyo_hist" +# f_in = sim + ".clm2.h1.2005-01.nc" +# f_out = sim + ".clm2.h1veg.0001-01.nc" +# f_dir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/DATA/" +# f_outfir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/new/OUTDIR/" +# +# cdf_file = nc.Dataset(f_dir + f_in, "r") +# +# ntim = cdf_file.variables['time'][:] +# nlat = cdf_file.variables['lat'][:] +# nlon = cdf_file.variables['lon'][:] +# +# grid1d_ixy = cdf_file.variables['grid1d_ixy'][:] +# grid1d_jxy = cdf_file.variables['grid1d_jxy'][:] +# grid1d_lon = cdf_file.variables['grid1d_lon'][:] +# grid1d_lat = cdf_file.variables['grid1d_lat'][:] +# land1d_lon = cdf_file.variables['land1d_lon'][:] +# land1d_lat = cdf_file.variables['land1d_lat'][:] +# land1d_ityplunit = cdf_file.variables['land1d_ityplunit'][:] +# pfts1d_lon = cdf_file.variables['pfts1d_lon'][:] +# pfts1d_lat = cdf_file.variables['pfts1d_lat'][:] +# pfts1d_active = cdf_file.variables['pfts1d_active'][:] +# pfts1d_itype_veg = cdf_file.variables['pfts1d_itype_veg'][:] +# pfts1d_wtgcell = cdf_file.variables['pfts1d_wtgcell'][:] +# pfts1d_wtlunit = cdf_file.variables['pfts1d_wtlunit'][:] +# +# GPP = cdf_file.variables['GPP'][:] +# +# cdf_file.close() +# +# out_file = nc.Dataset(f_outfir + f_out, "w") +# +# out_file.createDimension('time', None) +# out_file.createDimension('lat', len(nlat)) +# out_file.createDimension('lon', len(nlon)) +# gppGrass = out_file.createVariable( +# 'gppGrass', 'f4', ('time', 'lat', 'lon'), fill_value=1.e36) +# gppShrub = out_file.createVariable( +# 'gppShrub', 'f4', ('time', 'lat', 'lon'), fill_value=1.e36) +# gppTree = out_file.createVariable( +# 'gppTree', 'f4', ('time', 'lat', 'lon'), fill_value=1.e36) +# +# print 'Looking for grass' +# gppGrass[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'grass', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, +# grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, +# pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, +# pfts1d_wtgcell, pfts1d_wtlunit) +# print 'Looking for shrubs' +# gppShrub[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'shrub', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, +# grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, +# pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, +# pfts1d_wtgcell, pfts1d_wtlunit) +# print 'Looking for trees' +# gppTree[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'tree', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, +# grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, +# pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, +# pfts1d_wtgcell, pfts1d_wtlunit) +# +# out_file.close() +# +# +# if __name__ == '__main__': +# main() From 9ac793c20ee9d151562ecd58bb0cd61d04b7b0b5 Mon Sep 17 00:00:00 2001 From: sherimickelson Date: Mon, 16 Jul 2018 08:34:21 -0600 Subject: [PATCH 4/4] Modifications to comply with CMIP6 standards -Some minor modifications needed to comply with prepare -Revert back to older version of CLM_pft_toCMIP6_vegtype. It doesn't work with our default numpy. Will revert back later when we can verify we can run with newer version of numpy -Increase version number --- scripts/iconform | 30 +- source/pyconform/miptableparser.py | 1 + .../modules/CLM_pft_to_CMIP6_vegtype.py | 339 +++++++++++------- source/pyconform/modules/commonfunctions.py | 6 +- source/pyconform/version.py | 2 +- 5 files changed, 223 insertions(+), 155 deletions(-) diff --git a/scripts/iconform b/scripts/iconform index a712677a..70325a0f 100755 --- a/scripts/iconform +++ b/scripts/iconform @@ -166,10 +166,12 @@ def fill_missing_glob_attributes(attr, table, v, grids): attr["variable_id"] = v["variable_id"] if "branch_method" in attr.keys(): if "none" not in attr["branch_method"]: - #if "branch_time_in_child" in attr.keys(): - # attr["branch_time_in_child"] = "Get correct date format" - #if "branch_time_in_parent" in attr.keys(): - # attr["branch_time_in_parent"] = "Get correct date format" + if "branch_time_in_child" in attr.keys(): + if len(attr["branch_time_in_child"])>0: + attr["branch_time_in_child"] = float(attr["branch_time_in_child"].split('D')[0]) + if "branch_time_in_parent" in attr.keys(): + if len(attr["branch_time_in_parent"])>0: + attr["branch_time_in_parent"] = float(attr["branch_time_in_parent"].split('D')[0]) if "parent_mip_era" in attr.keys(): attr["parent_mip_era"] = attr["mip_era"] if "parent_source_id" in attr.keys(): @@ -190,13 +192,13 @@ def fill_missing_glob_attributes(attr, table, v, grids): if "variant_label" in attr.keys(): pre = attr["variant_label"].split('r')[1] - attr["realization_index"] = (pre.split('i')[0]) + attr["realization_index"] = int(pre.split('i')[0]) pre = pre.split('i')[1] - attr["initialization_index"] = (pre.split('p')[0]) + attr["initialization_index"] = int(pre.split('p')[0]) pre = pre.split('p')[1] - attr["physics_index"] = (pre.split('f')[0]) + attr["physics_index"] = int(pre.split('f')[0]) pre = int(pre.split('f')[1]) - attr["forcing_index"] = str(pre) + attr["forcing_index"] = int(pre) if "further_info_url" in attr.keys(): if "__FILL__" in attr["further_info_url"]: @@ -225,7 +227,7 @@ def fill_missing_glob_attributes(attr, table, v, grids): else: ripf = '' info_url = "{0}.{1}.{2}.{3}.{4}.{5}".format(mip_era, institution_id, source_id, experiment_id, sub_experiment_id, ripf) - attr['further_info_url'] = "http://furtherinfo.es-doc.org/" + info_url + attr['further_info_url'] = "https://furtherinfo.es-doc.org/" + info_url if "grid" in attr.keys(): if len(attr["realm"])>0: attr["grid"] = grids[attr["realm"].split()[0]] @@ -260,10 +262,10 @@ def defineVar(v, varName, attr, table_info, definition, ig, experiment, out_dir) # Get variables needed to piece together the filename ripf_list = ['realization_index','initialization_index','physics_index','forcing_index'] if all (ripf in attributes for ripf in ripf_list): - ripf = ("r{0}i{1}p{2}f{3}".format(attributes['realization_index'], - attributes['initialization_index'], - attributes['physics_index'], - attributes['forcing_index'])) + ripf = ("r{0}i{1}p{2}f{3}".format(str(attributes['realization_index']), + str(attributes['initialization_index']), + str(attributes['physics_index']), + str(attributes['forcing_index']))) else: ripf = '' @@ -515,7 +517,7 @@ def create_output(exp_dict, definitions, input_glob, attributes, output_path, ar with open(f, 'w') as outfile: json.dump(t, outfile, sort_keys=True, indent=4) else: - ignore = ['latitude','longitude','olevel','plev19','time','time1','alevhalf','ygre','xgre','vegtype','spectband','areacellg','alevel','xant','yant','rho','tau','plev3','gridlatitude','plev39','plev4','plev27','plev3','plev7h','plev7c','plev8','plev7h'] + ignore = ['latitude','longitude','olevel','plev19','time','time1','alevhalf','ygre','xgre','vegtype','spectband','areacellg','alevel','xant','yant','rho','tau','plev3','gridlatitude','plev39','plev4','plev27','plev3','plev7h','plev7c','plev8','plev7h','sdepth','siline','basin','olevel','site','soilpools','snowdepth','snowband','vegtype'] for n,t in TableSpec.iteritems(): for vn,var in t.iteritems(): if vn not in ignore: diff --git a/source/pyconform/miptableparser.py b/source/pyconform/miptableparser.py index 11359ef3..57df9745 100644 --- a/source/pyconform/miptableparser.py +++ b/source/pyconform/miptableparser.py @@ -421,6 +421,7 @@ def parse_table(self,exp,mips,tables,v_list,table_var_fields,table_axes_fields,t var['mipTable']=c_var.mipTable if c_var.mipTable in tables or '--ALL--' in tables: var["_FillValue"] = "1e+20" + #var["missing_value"] = "1e+20" if hasattr(c_var,'deflate'): var['deflate']= c_var.deflate if hasattr(c_var,'deflate_level'): diff --git a/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py b/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py index 759cc1ef..73875501 100644 --- a/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py +++ b/source/pyconform/modules/CLM_pft_to_CMIP6_vegtype.py @@ -1,8 +1,9 @@ #! /usr/bin/env python +import time, sys import numpy as np -from pyconform.physarray import PhysArray -from pyconform.functions import Function +from pyconform.physarray import PhysArray, UnitsError, DimensionsError +from pyconform.functions import Function, is_constant class CLM_pft_to_CMIP6_vegtype_Function(Function): @@ -10,17 +11,18 @@ class CLM_pft_to_CMIP6_vegtype_Function(Function): numargs = 18 def __init__(self, GPP, vegType, time, lat, lon, grid1d_ixy, grid1d_jxy, grid1d_lon, - grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, - pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, - pfts1d_wtgcell, pfts1d_wtlunit): + grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, + pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, + pfts1d_wtgcell, pfts1d_wtlunit): super(CLM_pft_to_CMIP6_vegtype_Function, self).__init__(GPP, vegType, time, lat, lon, grid1d_ixy, grid1d_jxy, grid1d_lon, - grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, - pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, - pfts1d_wtgcell, pfts1d_wtlunit) + grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, + pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, + pfts1d_wtgcell, pfts1d_wtlunit) def __getitem__(self, index): + pGPP = self.arguments[0][index] # vegType = grass, shrub, or tree vegType = self.arguments[1] @@ -30,22 +32,22 @@ def __getitem__(self, index): pgrid1d_ixy = self.arguments[5][index] pgrid1d_jxy = self.arguments[6][index] pgrid1d_lon = self.arguments[7][index] - pgrid1d_lat = self.arguments[8][index] + pgrid1d_lat = self.arguments[8][index] pland1d_lon = self.arguments[9][index] pland1d_lat = self.arguments[10][index] pland1d_ityplunit = self.arguments[11][index] - ppfts1d_lon = self.arguments[12][index] + ppfts1d_lon = self.arguments[12][index] ppfts1d_lat = self.arguments[13][index] ppfts1d_active = self.arguments[14][index] ppfts1d_itype_veg = self.arguments[15][index] - ppfts1d_wtgcell = self.arguments[16][index] + ppfts1d_wtgcell = self.arguments[16][index] ppfts1d_wtlunit = self.arguments[17][index] if index is None: - return PhysArray(np.zeros((0, 0, 0)), dimensions=[ptime.dimensions[0], plat.dimensions[0], plon.dimensions[0]]) - + return PhysArray(np.zeros((0,0,0)), dimensions=[ptime.dimensions[0],plat.dimensions[0],plon.dimensions[0]]) + GPP = pGPP.data - time = ptime.data + time = ptime.data lat = plat.data lon = plon.data grid1d_ixy = pgrid1d_ixy.data @@ -61,128 +63,191 @@ def __getitem__(self, index): pfts1d_itype_veg = ppfts1d_itype_veg.data pfts1d_wtgcell = ppfts1d_wtgcell.data pfts1d_wtlunit = ppfts1d_wtlunit.data - - # If 1, pft is active - active_pft = 1 - - # If 1, landunit is veg - veg_lunit = 1 - - # Veg-types and valid pfts - if 'tree' in vegType: - beg_pfts = 1 - end_pfts = 8 - elif 'shrub' in vegType: - beg_pfts = 9 - end_pfts = 11 - elif 'grass' in vegType: - beg_pfts = 12 - end_pfts = 14 - valid_pfts = ((pfts1d_wtgcell > 0.) * - (pfts1d_itype_veg >= beg_pfts) * - (pfts1d_itype_veg <= end_pfts)) - - # Will contain weighted average for grass pfts on 2d grid - varo_vegType = np.ones([len(time), len(lat), len(lon)]) * 1e20 - - vegtyp_lunit_indices = np.where(land1d_ityplunit == veg_lunit) - active_ptype_indices = np.where(pfts1d_active == active_pft) - - gcell_pts = np.stack((grid1d_lon, grid1d_lat), axis=1) - lunit_pts = np.stack( - (land1d_lon[vegtyp_lunit_indices], land1d_lat[vegtyp_lunit_indices]), axis=1) - pfts_pts = np.stack( - (pfts1d_lon[active_ptype_indices], pfts1d_lat[active_ptype_indices]), axis=1) - - gcell_indices = np.where(np.all(np.isin(gcell_pts, lunit_pts), axis=1)) - gcell_pts = gcell_pts[gcell_indices] - gcell_indices = np.where(np.all(np.isin(gcell_pts, pfts_pts), axis=1)) - - gcell_lon = grid1d_lon[gcell_indices] - gcell_lat = grid1d_lat[gcell_indices] - gcell_ixy = grid1d_ixy[gcell_indices] - gcell_jxy = grid1d_jxy[gcell_indices] - - for lon, lat, i, j in np.nditer([gcell_lon, gcell_lat, gcell_ixy, gcell_jxy]): - pft_idx = np.where( - np.all((lon, lat) == pfts_pts, axis=1) * valid_pfts)[0] - dum = GPP[:, pft_idx] - pfts1d_wtlunit_veg = (pfts1d_wtlunit[pft_idx]).astype(np.float32) - weights = pfts1d_wtlunit_veg / np.sum(pfts1d_wtlunit_veg) - varo_vegType[:, j - 1, i - 1] = np.sum(dum * weights) - - new_name = 'CLM_pft_to_CMIP6_vegtype({},{})'.format(pGPP.name, vegType) - - varo_vegType[varo_vegType >= 1e+16] = 1e+20 + + + # vegType = grass, shrub, or tree + + # Tolerance check for weights summing to 1 + eps = 1.e-5 + + # If 1, pft is active + active_pft = 1 + + # If 1, landunit is veg + veg_lunit = 1 + + # C3 arctic grass, + # C3 non-arctic grass, + # C4 grass + beg_grass_pfts = 12 + end_grass_pfts = 14 + + # broadleaf evergreen shrub - temperate, + # broadleaf deciduous shrub - temperate, + # broadleaf deciduous shrub - boreal + beg_shrub_pfts = 9 + end_shrub_pfts = 11 + + # needleleaf evergreen tree - temperate, + # needleleaf evergreen tree - boreal, + # needleleaf deciduous tree - boreal, + # broadleaf evergreen tree - tropical, + # broadleaf evergreen tree - temperate, + # broadleaf deciduous tree - tropical, + # broadleaf deciduous tree - temperate, + # broadleaf deciduous tree - boreal + beg_tree_pfts = 1 + end_tree_pfts = 8 + + # Will contain weighted average for grass pfts on 2d grid + varo_vegType = np.zeros([len(time),len(lat),len(lon)]) + tu = np.stack((pfts1d_lon,pfts1d_lat, pfts1d_active), axis=1) + + ind = np.stack((grid1d_ixy,grid1d_jxy), axis=1) + + lu = np.stack((land1d_lon,land1d_lat,land1d_ityplunit), axis=1) + + # Loop over lat/lons + for ixy in range(len(lon)): + for jxy in range(len(lat)): + grid_indx = -99 + # 1d grid index + ind_comp = (ixy+1,jxy+1) + gi = np.where(np.all(ind==ind_comp, axis=1))[0] + if len(gi) > 0: + grid_indx = gi[0] + + # Check for valid land gridcell + if grid_indx != -99: + + # Gridcell lat/lons + grid1d_lon_pt = grid1d_lon[grid_indx] + grid1d_lat_pt = grid1d_lat[grid_indx] + + # veg landunit index for this gridcell + t_var = (grid1d_lon_pt, grid1d_lat_pt, veg_lunit) + landunit_indx = np.where(np.all(t_var == lu, axis=1))[0] + + # Check for valid veg landunit + if landunit_indx.size > 0: + if 'grass' in vegType: + t_var = (grid1d_lon_pt,grid1d_lat_pt, active_pft) + pft_indx = np.where( np.all(t_var == tu, axis=1) * (pfts1d_wtgcell > 0.) * (pfts1d_itype_veg >= beg_grass_pfts) * (pfts1d_itype_veg <= end_grass_pfts))[0] + elif 'shrub' in vegType: + t_var = (grid1d_lon_pt,grid1d_lat_pt, active_pft) + pft_indx = np.where( np.all(t_var==tu, axis=1) * (pfts1d_wtgcell > 0.) * (pfts1d_itype_veg >= beg_shrub_pfts) * (pfts1d_itype_veg <= end_shrub_pfts))[0] + elif 'tree' in vegType: + t_var = (grid1d_lon_pt,grid1d_lat_pt, active_pft) + pft_indx = np.where( np.all(t_var==tu, axis=1) * (pfts1d_wtgcell > 0.) * (pfts1d_itype_veg >= beg_tree_pfts) * (pfts1d_itype_veg <= end_tree_pfts))[0] + + # Check for valid pfts and compute weighted average + if pft_indx.size > 0: + if 'grass' in vegType: + pfts1d_wtlunit_grass = (pfts1d_wtlunit[pft_indx]).astype(np.float32) + dum = GPP[:,pft_indx] + weights = pfts1d_wtlunit_grass / np.sum(pfts1d_wtlunit_grass) + if np.absolute(1.-np.sum(weights)) > eps: + print("Weights do not sum to 1, exiting") + sys.exit(-1) + varo_vegType[:,jxy,ixy] = np.sum(dum * weights) + + elif 'shrub' in vegType: + pfts1d_wtlunit_shrub = (pfts1d_wtlunit[pft_indx]).astype(np.float32) + dum = GPP[:,pft_indx] + weights = pfts1d_wtlunit_shrub / np.sum(pfts1d_wtlunit_shrub) + varo_vegType[:,jxy,ixy] = np.sum(dum * weights) + + elif 'tree' in vegType: + pfts1d_wtlunit_tree = (pfts1d_wtlunit[pft_indx]).astype(np.float32) + dum = GPP[:,pft_indx] + weights = pfts1d_wtlunit_tree / np.sum(pfts1d_wtlunit_tree) + varo_vegType[:,jxy,ixy] = np.sum(dum * weights) + + else: + varo_vegType[:,jxy,ixy] = 1e+20 + else: + varo_vegType[:,jxy,ixy] = 1e+20 + else: + varo_vegType[:,jxy,ixy] = 1e+20 + + + new_name = 'CLM_pft_to_CMIP6_vegtype({}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{})'.format(pGPP.name, + vegType, ptime.name, plat.name, plon.name, pgrid1d_ixy.name, pgrid1d_jxy.name, pgrid1d_lon.name, + pgrid1d_lat.name, pland1d_lon.name, pland1d_lat.name, pland1d_ityplunit.name, + ppfts1d_lon.name, ppfts1d_lat.name, ppfts1d_active.name, ppfts1d_itype_veg.name, + ppfts1d_wtgcell.name, ppfts1d_wtlunit.name) + + print 'FINISHED FUNCTION' + + varo_vegType[varo_vegType>=1e+16] = 1e+20 ma_varo_vegType = np.ma.masked_values(varo_vegType, 1e+20) - return PhysArray(ma_varo_vegType, name=new_name, units=pGPP.units) - - -# def main(): -# -# import netCDF4 as nc -# -# sim = "clm50_r243_1deg_GSWP3V2_cropopt_nsc_emergeV2F_dailyo_hist" -# f_in = sim + ".clm2.h1.2005-01.nc" -# f_out = sim + ".clm2.h1veg.0001-01.nc" -# f_dir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/DATA/" -# f_outfir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/new/OUTDIR/" -# -# cdf_file = nc.Dataset(f_dir + f_in, "r") -# -# ntim = cdf_file.variables['time'][:] -# nlat = cdf_file.variables['lat'][:] -# nlon = cdf_file.variables['lon'][:] -# -# grid1d_ixy = cdf_file.variables['grid1d_ixy'][:] -# grid1d_jxy = cdf_file.variables['grid1d_jxy'][:] -# grid1d_lon = cdf_file.variables['grid1d_lon'][:] -# grid1d_lat = cdf_file.variables['grid1d_lat'][:] -# land1d_lon = cdf_file.variables['land1d_lon'][:] -# land1d_lat = cdf_file.variables['land1d_lat'][:] -# land1d_ityplunit = cdf_file.variables['land1d_ityplunit'][:] -# pfts1d_lon = cdf_file.variables['pfts1d_lon'][:] -# pfts1d_lat = cdf_file.variables['pfts1d_lat'][:] -# pfts1d_active = cdf_file.variables['pfts1d_active'][:] -# pfts1d_itype_veg = cdf_file.variables['pfts1d_itype_veg'][:] -# pfts1d_wtgcell = cdf_file.variables['pfts1d_wtgcell'][:] -# pfts1d_wtlunit = cdf_file.variables['pfts1d_wtlunit'][:] -# -# GPP = cdf_file.variables['GPP'][:] -# -# cdf_file.close() -# -# out_file = nc.Dataset(f_outfir + f_out, "w") -# -# out_file.createDimension('time', None) -# out_file.createDimension('lat', len(nlat)) -# out_file.createDimension('lon', len(nlon)) -# gppGrass = out_file.createVariable( -# 'gppGrass', 'f4', ('time', 'lat', 'lon'), fill_value=1.e36) -# gppShrub = out_file.createVariable( -# 'gppShrub', 'f4', ('time', 'lat', 'lon'), fill_value=1.e36) -# gppTree = out_file.createVariable( -# 'gppTree', 'f4', ('time', 'lat', 'lon'), fill_value=1.e36) -# -# print 'Looking for grass' -# gppGrass[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'grass', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, -# grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, -# pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, -# pfts1d_wtgcell, pfts1d_wtlunit) -# print 'Looking for shrubs' -# gppShrub[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'shrub', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, -# grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, -# pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, -# pfts1d_wtgcell, pfts1d_wtlunit) -# print 'Looking for trees' -# gppTree[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'tree', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, -# grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, -# pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, -# pfts1d_wtgcell, pfts1d_wtlunit) -# -# out_file.close() -# -# -# if __name__ == '__main__': -# main() + return PhysArray(ma_varo_vegType, name=new_name, units=pGPP.units) + + +def main(argv=None): + + import netCDF4 as nc + + sim = "clm50_r243_1deg_GSWP3V2_cropopt_nsc_emergeV2F_dailyo_hist" + f_in = sim+".clm2.h1.2005-01.nc" + f_out = sim+".clm2.h1veg.0001-01.nc" + f_dir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/DATA/" + f_outfir = "/glade2/scratch2/mickelso/CMIP6_LND_SCRIPTS/new/OUTDIR/" + + cdf_file = nc.Dataset(f_dir+f_in,"r") + + ntim = cdf_file.variables['time'][:] + nlat = cdf_file.variables['lat'][:] + nlon = cdf_file.variables['lon'][:] + + grid1d_ixy = cdf_file.variables['grid1d_ixy'][:] + grid1d_jxy = cdf_file.variables['grid1d_jxy'][:] + grid1d_lon = cdf_file.variables['grid1d_lon'][:] + grid1d_lat = cdf_file.variables['grid1d_lat'][:] + land1d_lon = cdf_file.variables['land1d_lon'][:] + land1d_lat = cdf_file.variables['land1d_lat'][:] + land1d_ityplunit = cdf_file.variables['land1d_ityplunit'][:] + pfts1d_lon = cdf_file.variables['pfts1d_lon'][:] + pfts1d_lat = cdf_file.variables['pfts1d_lat'][:] + pfts1d_active = cdf_file.variables['pfts1d_active'][:] + pfts1d_itype_veg = cdf_file.variables['pfts1d_itype_veg'][:] + pfts1d_wtgcell = cdf_file.variables['pfts1d_wtgcell'][:] + pfts1d_wtlunit = cdf_file.variables['pfts1d_wtlunit'][:] + + + GPP = cdf_file.variables['GPP'][:] + + cdf_file.close() + + out_file = nc.Dataset(f_outfir+f_out,"w") + + time = out_file.createDimension('time',None) + lat = out_file.createDimension('lat',len(nlat)) + lon = out_file.createDimension('lon',len(nlon)) + gppGrass = out_file.createVariable('gppGrass', 'f4', ('time', 'lat', 'lon'),fill_value=1.e36) + gppShrub = out_file.createVariable('gppShrub', 'f4', ('time', 'lat', 'lon'),fill_value=1.e36) + gppTree = out_file.createVariable('gppTree', 'f4', ('time', 'lat', 'lon'),fill_value=1.e36) + + print 'Looking for grass' + gppGrass[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'grass', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, + grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, + pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, + pfts1d_wtgcell, pfts1d_wtlunit) + print 'Looking for shrubs' + gppShrub[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'shrub', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, + grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, + pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, + pfts1d_wtgcell, pfts1d_wtlunit) + print 'Looking for trees' + gppTree[:] = CLM_pft_to_CMIP6_vegtype(GPP, 'tree', ntim, nlat, nlon, grid1d_ixy, grid1d_jxy, grid1d_lon, + grid1d_lat, land1d_lon, land1d_lat, land1d_ityplunit, + pfts1d_lon, pfts1d_lat, pfts1d_active, pfts1d_itype_veg, + pfts1d_wtgcell, pfts1d_wtlunit) + + + out_file.close() + +if __name__ == '__main__': + main() + diff --git a/source/pyconform/modules/commonfunctions.py b/source/pyconform/modules/commonfunctions.py index b9501b05..51693e85 100644 --- a/source/pyconform/modules/commonfunctions.py +++ b/source/pyconform/modules/commonfunctions.py @@ -525,11 +525,11 @@ def __getitem__(self, index): #=================================================================================================== # soilpoolsFunction #=================================================================================================== -class soilpoolsFunction(Function): - key = 'soilpools' +class get_soilpoolsFunction(Function): + key = 'get_soilpools' def __init__(self, p_data1,p_data2,p_data3,p_soilpool): - super(soilpoolsFunction, self).__init__(p_data1,p_data2,p_data3,p_soilpool) + super(get_soilpoolsFunction, self).__init__(p_data1,p_data2,p_data3,p_soilpool) def __getitem__(self, index): p_data1 = self.arguments[0][index] diff --git a/source/pyconform/version.py b/source/pyconform/version.py index ce5f98f0..d0c9d782 100644 --- a/source/pyconform/version.py +++ b/source/pyconform/version.py @@ -1,2 +1,2 @@ # Single place for version information -__version__ = '0.2.5' +__version__ = '0.2.6'