農作物

def tif2nc(tif_name,nc_name,lev):
  import numpy as np
  import netCDF4
  from pyproj import Proj
  import rasterio
  import numpy as np

  img = rasterio.open(tif_name)
  nx,ny,nz=img.width,img.height,img.count
  if nz!=1:return -1
  dx,dy=360./(nx-1),180./(ny-1)
  lon_1d=[-180+dx*i for i in range(nx)]
  lat_1d=[90-dy*i for i in range(ny)]
  data=img.read()
  lonm, latm = np.meshgrid(lon_1d, lat_1d)
  DD={'lon':lonm.flatten(),'lat':latm.flatten(),'val':data.flatten()}
  df=DataFrame(DD)
  boo=(df.lon>=60)&(df.lon<=180)&(df.lat>=-10)&(df.lat<=50)
  df=df.loc[boo].reset_index(drop=True)
  
  nc = netCDF4.Dataset(nc_name, 'r+')
  pnyc = Proj(proj='lcc', datum='NAD83', lat_1=nc.P_ALP, lat_2=nc.P_BET,lat_0=nc.YCENT, lon_0=nc.XCENT, x_0=0, y_0=0.0)
  V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
  nt,nlay,nrow,ncol=(nc.variables[V[3][0]].shape[i] for i in range(4))

  #d00範圍：北緯-10~50、東經60~180。'area': [50, 60, -10, 180,],
  x,y=pnyc(list(df.lon),list(df.lat), inverse=False)
  x,y=np.array(x),np.array(y)
  df['ix']=np.array((x-nc.XORIG)/nc.XCELL,dtype=int)
  df['iy']=np.array((y-nc.YORIG)/nc.YCELL,dtype=int)
  df=df.loc[(df.ix>=0)&(df.ix<ncol)&(df.iy>=0)&(df.iy<nrow)].reset_index(drop=True)
  df['ixy']=[str(i)+'_'+str(j) for i,j in zip(df.ix,df.iy)]
  pv=pivot_table(df,index='ixy',values='val',aggfunc=np.sum).reset_index()
  pv['ix']=[int(i.split('_')[0]) for i in pv.ixy]
  pv['iy']=[int(i.split('_')[1]) for i in pv.ixy]
  var=np.zeros(shape=(nrow,ncol))
  var[pv.iy,pv.ix]=pv.val
  if lev==3: #beansedible
    nc[V[3][0]][0,lev,:,:]+=var[:,:]
  else:
    nc[V[3][0]][0,lev,:,:]=0.
    if lev==19: #wheat_spring
      nc[V[3][0]][0,lev,320:,:]=var[320:,:]
    elif lev==20: #wheat_winter
      nc[V[3][0]][0,lev,:320,:]=var[:320,:]
    else:
      nc[V[3][0]][0,lev,:,:]+=var[:,:]
  nc.close()
  return 0

nc_name='temp.nc'
fnameO='d21_175.json'
with open(fnameO,'r') as jsonfile:
  d21_175=json.load(jsonfile)
crop21=list(d21_175)
crop21.sort()
#18項1對1
s=set([d21_175[i] for i in crop21 if i not in ['wheat_spring','wheat_winter','beansedible']]) 
d175_21={d21_175[i]:crop21.index(i) for i in s}
fn1,fn2,fn3='HarvestedAreaYield175Crops_Geotiff/','_HarvAreaYield_Geotiff/','_HarvestedAreaHectares.tif'
for c in s:
  tif_name=fn1+c+fn2+c+fn3
  i=tif2nc(tif_name,nc_name,d175_21[c])
#crop21.index('beansedible')=3，多對1
for c in ['broadbean', 'greenbean', 'greenbroadbean', 'stringbean']:
  tif_name=fn1+c+fn2+c+fn3
  i=tif2nc(tif_name,nc_name,3)
#1對多
c='wheat' 
tif_name=fn1+c+fn2+c+fn3
for lev in [19,20]:
  i=tif2nc(tif_name,nc_name,lev)

fname='irr28.nc'
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]]
irr=nc[V[3][0]][0,0,:,:] #28種灌溉形式之總合(area fraction)

fname='temp.nc'
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]]
v=V[3][0]
nt,nlay,nrow,ncol=(nc.variables[v].shape[i] for i in range(4))
km2=nc[v][0,:21,:,:]*0.01
for lev in range(21,42):
  nc[v][0,lev,:,:]=irr[:,:]*km2[lev-21,:,:]
for lev in range(21):
  nc[v][0,lev,:,:]=(1-irr[:,:])*km2[lev,:,:]
svar=np.sum(nc[v][0,:,:,:],axis=0)
a=np.where(svar<255,255,svar)
nc[v][0,:,:,:]=nc[v][0,:,:,:]/a[None,:,:]
nc.close()