衛星數據轉換成D4座標系統

Table of contents

背景
程式說明

背景

此處將Sinusoidal Tile Grid檔案之經緯度系統，轉成CMAQ之D4範圍座標系統。
因二者的解析度相同，沒有必要重新進行內插或regrid計算，直接由陣列中選取最接近的點，移轉其數值到新的系統即可，作法如建立cmaq邊界點位置與cams網格系統的對照關係或LGHAP數據之切割與平均。
麻煩之處是台灣正好處在Sinusoidal Tile Grid2格之間，因此需從2格共2,880,000個點當中，找到符合D4範圍的座標點，建立其對應關係，以避免矩陣太大難以處理。

程式說明

IO’s

引數：無
'2016allnc'+hh+'.txt'，hh=['28','29']，將從其中各讀取1個檔案作為Sinusoidal Tile Grid格點經緯度之來源。
tempTW.nc：為CMAQ D4範圍(解析度1公里)之模板，內容如下

netcdf tempTW {
dimensions:
        TSTEP = UNLIMITED ; // (24 currently)
        VAR = 1 ;
        DATE-TIME = 2 ;
        LAY = 4 ;
        ROW = 393 ;
        COL = 276 ;
variables:
        int TFLAG(TSTEP, VAR, DATE-TIME) ;
...
                :GDTYP = 2 ;
                :P_ALP = 10. ;
                :P_BET = 40. ;
                :P_GAM = 120. ;
                :XCENT = 120. ;
                :YCENT = 25. ;
                :XORIG = -72000. ;
                :YORIG = -345000. ;
                :NCOLS = 276L ;
                :NROWS = 393L ;
                :XCELL = 1000L ;
                :YCELL = 1000L ;
                :NLAYS = 4L ;
                :VGLVLS = 1., 0.995, 0.5, 0.2, 0.1 ;

N393_276.bin：結果檔案，將記存者新座標(D4)對應到舊座標(V06H28~29)之引數。

新座標系統之計算

D4為環保署公版模示範圍，其網格已經是正交直角系統，只需簡單計算即可。
按照新座標(D4)定義pyproj.Proj備用。

# open a blank template nc file for D4 domain
fname='/nas2/cmaqruns/2019TZPP/output/Annual/aTZPP/MCD19A2.006/tempTW.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]]
nt,nlay,nrow,ncol=(nc[V[3][0]].shape[i] for i in range(4))

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)
X=[nc.XORIG+nc.XCELL*i for i in range(ncol)]
Y=[nc.YORIG+nc.YCELL*i for i in range(nrow)]
x,y=np.meshgrid(X, Y)
x=x.flatten()
y=y.flatten()
nbnd1=len(x)
minx=x[0]-5000;maxx=x[-1]+5000;miny=y[0]-5000;maxy=y[-1]+5000

舊座標系統(Sinusoidal Tile Grid)之計算

以np.where進行第1次篩選，選出D4範圍(含外圈5公里)之衛星數據格點。
- 因2個網格的點數共有1200 ×1200 ×2近3百萬個1公里格點，處理非常困難，因此需先進行篩檢。

lat=[];lon=[]
for hh in ['28','29']:
  with open('2016allnc'+hh+'.txt','r') as f:
    fnames=[i.strip('\n') for i in f]
  nc = netCDF4.Dataset(fnames[0],'r')
  lat=lat+list(np.array(nc['grid1km_latitude'][:,:]).flatten())
  lon=lon+list(np.array(nc['grid1km_longitude'][:,:]).flatten())
x1,y1=pnyc(lon,lat, inverse=False)
x1=np.array(x1)
y1=np.array(y1)
boo=(x1>=minx)&(y1>=miny)&(x1<=maxx)&(y1<=maxy)
idx0=np.where(boo)
x10=x1[idx0[0]]
y10=y1[idx0[0]]

找到並記錄最近點

每個新格點迴圈
每個新格點進行篩選，只計算附近5公里範圍的舊格點(第2次篩選)
找到其中最近點的位置(第3次篩選)
回復到原始的引數需進行3次套疊(n[i]=idx0[0][idx1[0][idx]])

n=[-1 for i in range(nbnd1)]
for i in range(nbnd1):
  minx=x[i]-5000;maxx=x[i]+5000
  miny=y[i]-5000;maxy=y[i]+5000
  boo=(x10>=minx)&(y10>=miny)&(x10<=maxx)&(y10<=maxy)
  idx1=np.where(boo)
  if len(idx1[0])==0:continue
  x11=x10[idx1[0]]
  y11=y10[idx1[0]]
  dist=(x[i]-x11)**2+(y[i]-y11)**2      #nearest grib data for bcon
  idx=np.where(dist==np.min(dist))[0][0]
  if type(idx)==list and len(idx)>1: idx=idx[0]
  n[i]=idx0[0][idx1[0][idx]]
nar=np.array(n,dtype=int)
fnameO='N393_276.bin'
with FortranFile(fnameO, 'w') as f:
    f.write_record(nar)

程式碼下載

Download: 切割並轉檔 genN_D4T.py