亞洲土地使用檔案

fname='beld4.CWBWRF_15k.ncf'
nc = netCDF4.Dataset(fname, 'r+')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
nrow,nt,nlay,ncol=(nc.variables[V[3][0]].shape[i] for i in range(4)
for j in range(nrow,389):
  for v in V[3]:
    nc[v][j,0,0,:]=0.
nc.NROWS=389
nc.close()

nc = netCDF4.Dataset(fname, 'r+')
for j in range(ncol,665):
  for v in V[3]:
    nc[v][j,0,0,:]=0.
nc.NCOLS=665
# 其他全域屬性，抄自mcip結果
# 注意：VG*、NLAYS等垂直設定不能抄襲、NVARS也不能覆蓋、刪除舊的HISTORY會使檔案看起來較為清爽
fname='../mcip/1804_run5/CWBWRF_15k/LUFRAC_CRO_1804_run5.nc'
nc0 = netCDF4.Dataset(fname, 'r')
atts=['CDATE',  'CTIME', 'EXEC_ID', 'FILEDESC', 'FTYPE', 'GDNAM', 'GDTYP', 'HISTORY', 'IOAPI_VERSION', 'NCO', 'NCOLS',  'NROWS',
     'NTHIK', 'P_ALP', 'P_BET', 'P_GAM', 'UPNAM', 'WDATE',
     'WTIME', 'XCELL', 'XCENT', 'XORIG', 'YCELL', 'YCENT', 'YORIG']
for i in atts:
  if i not in dir(nc0):continue
  exec('nc.'+i+'=nc0.'+i)
nc.close()

fname='beld4.CWBWRF_15k.ncf'
nc = netCDF4.Dataset(fname, 'r+')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
modis=["MODIS_"+str((i+1)%17) for i in range(17)]
for v in V[3]:
  nc[v][:]=0.
  if v[:6]=="MODIS_":
    try:
      ii=int(v[6:8])-1
    except:
      ii=17+int(v[9:10])
    else:
      if ii<0:ii=ii+17
    nc[v][0,0,:,:]=nc0["LUFRAC"][0,ii,:,:]

modis=[i for i in V[3] if 'MODIS' in i and 'Res' not in i]
var=np.zeros(shape=(len(modis),nrow,ncol))
for v in modis:
  iv=modis.index(v)
  var[iv,:,:]=nc[v][0,0,:,:]
sumv=np.sum(var[:,:,:],axis=0)
np.max(sumv) 
#1.0000000596046448