python轉換wrf輸出的資料為網頁視覺化json格式

2022-09-28 14:00:50

前言

一般網頁視覺化風場中的資料都是json格式，而如果我們希望將wrf模式模擬輸出的風場資料在網頁中進行展示，這就需要先將wrfoutput資料轉換為網頁可以識別的json格式。
這裡主要需要用到json庫，主要的實現方式就是將讀取的風場風量U，V轉換為字典並存到json檔案中
同時，由於wrf模擬的資料一般是非等間距的網格，需要先將資料進行插值，插值到等間距的網格，這裡可以通過NCL的函數rcm2rgrid_Wrap實現

舉個例子，將模式中設定為蘭伯特投影的網格：

插值為等間距網格：

主要的程式設計分為兩部分：

第一部分通過NCL指令碼將wrfout資料轉換為等間距網格，並匯出為netcdf格式；
第二部分通過python指令碼將第一步匯出的nc格式進行轉換，並儲存輸出為json格式。

NCL插值指令碼1

需要修改的就是路徑和變數，我下面展示指令碼不僅有風場資料u，v還有降水，海表面壓力，氣溫等，可自行修改

begin
  a = addfile("/Users/WRF/outdata/2022071000/wrfout_d01_2022-07-10_01:00:00","r")
  lat2d = a->XLAT(0,:,:)
  lon2d = a->XLONG(0,:,:)
  lat1d = lat2d(:,0)
  lon1d = lon2d(0,:)
 
  time = wrf_user_getvar(a,"XTIME",-1)
  u10 = wrf_user_getvar(a,"U10",0)
  v10 = wrf_user_getvar(a,"V10",0)
  slp = wrf_user_getvar(a,"slp",0)
  t2  = wrf_user_getvar(a,"T2",0)
  td  = wrf_user_getvar(a,"td",0)
  rainc = wrf_user_getvar(a,"RAINC",0)
  rainnc = wrf_user_getvar(a,"RAINNC",0)
  
  u10@lat2d = lat2d
  u10@lon2d = lon2d
  u10_ip = rcm2rgrid_Wrap(lat2d,lon2d,u10,lat1d,lon1d,0)
  
  v10@lat2d = lat2d
  v10@lon2d = lon2d
  v10_ip = rcm2rgrid_Wrap(lat2d,lon2d,v10,lat1d,lon1d,0)
  
  slp_ip  =   rcm2rgrid_Wrap(lat2d,lon2d,slp,lat1d,lon1d,0)
  t2_ip  =   rcm2rgrid_Wrap(lat2d,lon2d,t2,lat1d,lon1d,0)
  td_ip  =   rcm2rgrid_Wrap(lat2d,lon2d,td,lat1d,lon1d,0)
  rainc_ip  =   rcm2rgrid_Wrap(lat2d,lon2d,rainc,lat1d,lon1d,0)
  rainnc_ip  =   rcm2rgrid_Wrap(lat2d,lon2d,rainnc,lat1d,lon1d,0)
  outf = addfile("/Users/wrfout_d01_2022-07-10_01:00:00.nc","c")
  outf->time =  time
  outf->lat  =  lat2d
  outf->lon  =  lon2d
  outf->u10  =  u10_ip
  outf->v10  =  v10_ip
  outf->slp  =  slp_ip
  outf->t2   =  t2_ip
  outf->td   =  td_ip
  outf->rainc   =  rainc_ip
  outf->rainnc  =  rainnc_ip
end

上述指令碼的缺點在於只能基於模式模擬的經緯度區域進行插值，意思就是說他的經緯度區域是固定的那麼大

NCL插值指令碼2

NCL還有一個函數可以實現上述過程，就是ESMF_regrid，該函數的優點在於可以實現任意經緯度範圍的插值，但是不足在於對於存在高度層的變數，暫時無法進行高度層的資料讀取。

（也可能我水平有限不知道。。。。）這裡也附上指令碼：

load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

begin
  a = addfile("/Users/WRF/outdata/2022071000/wrfout_d01_2022-07-10_01:00:00","r")
  u10 = wrf_user_getvar(a,"U10",0)
  v10 = wrf_user_getvar(a,"V10",0)
  slp = wrf_user_getvar(a,"slp",0)
  t2  = wrf_user_getvar(a,"T2",0)
;  td  = wrf_user_getvar(a,"td",0)
  rainc = wrf_user_getvar(a,"RAINC",0)
  rainnc = wrf_user_getvar(a,"RAINNC",0)
  
  u10@lat2d = a->XLAT(0,:,:) 
  u10@lon2d = a->XLONG(0,:,:)
  v10@lat2d = a->XLAT(0,:,:) 
  v10@lon2d = a->XLONG(0,:,:)
  slp@lat2d = a->XLAT(0,:,:) 
  slp@lon2d = a->XLONG(0,:,:)
  t2@lat2d = a->XLAT(0,:,:) 
  t2@lon2d = a->XLONG(0,:,:)
;  td@lat2d = a->XLAT(0,:,:) 
;  td@lon2d = a->XLONG(0,:,:)
  rainc@lat2d = a->XLAT(0,:,:) 
  rainc@lon2d = a->XLONG(0,:,:)
  rainnc@lat2d = a->XLAT(0,:,:) 
  rainnc@lon2d = a->XLONG(0,:,:)
  
  lat2d = a->XLAT(0,:,:)
  lon2d = a->XLONG(0,:,:)
  lat1d = lat2d(:,0)
  lon1d = lon2d(0,:)
  latS = -20
  latN = 50
  lonW = 95
  lonE = 145

  Opt = True
  Opt@InterpMethod = "bilinear" 
  Opt@ForceOverwrite = True 
  
  Opt@SrcMask2D = where(.not. ismissing(v10),1,0) 
  Opt@DstGridType = "0.1deg"
  Opt@DstLLCorner = (/latS, lonW /) 
  Opt@DstURCorner = (/latN, lonE /) 
  
  u10_regrid = ESMF_regrid(u10,Opt)
  v10_regrid = ESMF_regrid(v10,Opt)
  slp_regrid = ESMF_regrid(slp,Opt)
  t2_regrid = ESMF_regrid(t2,Opt)
;  td_regrid = ESMF_regrid(td,Opt)
  rainc_regrid = ESMF_regrid(rainc,Opt)
  rainnc_regrid = ESMF_regrid(rainnc,Opt)
  
  time = wrf_user_getvar(a,"XTIME",-1)
  
  nlon = dimsizes(v10_regrid&lon)
  nlat = dimsizes(v10_regrid&lat)
  
  ofile = "wrfout_d01_2022-07-10_01:00:00.nc"
  system("rm -rf "+ofile) 
  fout = addfile(ofile,"c") 
  
  dimNames = (/"lat", "lon"/)
  dimSizes = (/nlat, nlon/)
  dimUnlim = (/False, False/)
  
  filedimdef(fout,dimNames,dimSizes,dimUnlim) ;-- define dimensions
  
  filevardef(fout,"lat",typeof(v10_regrid&lat),getvardims(v10_regrid&lat))
  filevardef(fout,"lon",typeof(v10_regrid&lon),getvardims(v10_regrid&lon))
  
  filevardef(fout,"u10",typeof(u10_regrid),getvardims(u10_regrid))
  filevardef(fout,"v10",typeof(v10_regrid),getvardims(v10_regrid))
  filevardef(fout,"slp",typeof(slp_regrid),getvardims(slp_regrid))
  filevardef(fout,"t2",typeof(t2_regrid),getvardims(t2_regrid))
;  filevardef(fout,"td",typeof(td_regrid),getvardims(td_regrid))
  filevardef(fout,"rainc",typeof(rainc_regrid),getvardims(rainc_regrid))
  filevardef(fout,"rainnc",typeof(rainnc_regrid),getvardims(rainnc_regrid))
    
  filevarattdef(fout,"lat",v10_regrid&lat) ;-- copy lat attributes
  filevarattdef(fout,"lon",v10_regrid&lon) ;-- copy lon attributes
  filevarattdef(fout,"u10",u10_regrid)
  filevarattdef(fout,"v10",v10_regrid)
  filevarattdef(fout,"slp",slp_regrid)
  filevarattdef(fout,"t2",t2_regrid)
;  filevarattdef(fout,"td",td_regrid)
  filevarattdef(fout,"rainc",rainc_regrid)
  filevarattdef(fout,"rainnc",rainnc_regrid)
  
  setfileoption(fout,"DefineMode",False)

  fout->u10 = (/u10_regrid/)
  fout->v10 = (/v10_regrid/) 
  fout->slp = (/slp_regrid/) 
  fout->t2 = (/t2_regrid/) 
;  fout->td = (/td_regrid/) 
  fout->rainc  = (/rainc_regrid/) 
  fout->rainnc = (/rainnc_regrid/) 
  
  fout->lat = (/v10_regrid&lat/) ;-- write lat to new netCDF file
  fout->lon = (/v10_regrid&lon/) ;-- write lon to new netCDF file
  fout->time =  time
end

PS:執行該指令碼會生成四個nc檔案，分別為：destination_grid_file.nc、source_grid_file.nc、weights_file.nc、wrfout_d01_2022-07-10_01:00:00.nc。其中，wrfout_d01_2022-07-10_01:00:00.nc是我需要的檔案，但是其他三個檔案如何在執行指令碼的過程去掉暫未解決。

python格式轉換指令碼1

python指令碼如下所示：

# -*- coding: utf-8 -*-
"""
Created on %(date)s

@author: %(jixianpu)s

Email : 211311040008@hhu.edu.cn

introduction : keep learning althongh walk slowly
"""
"""
用來讀取用ncl插值後的wrfoutput.nc 資料,並生成對應檔名的json格式
"""
import pandas as pd
import os
import json
import netCDF4 as nc
import numpy as np
import  datetime
from netCDF4 import Dataset
import argparse
from argparse import RawDescriptionHelpFormatter
import xarray as xr
import sys
import glob

date = sys.argv[1]
date = str(date)
frst = sys.argv[2]
step = sys.argv[3]

path = r'/Users/WRF/outdata/2022071000/'#只能是已經存在的檔案目錄且有資料才可以進行讀取
start = datetime.datetime.strptime(date,'%Y%m%d%H').strftime("%Y-%m-%d_%H:%M:%S")
end = (datetime.datetime.strptime(date,'%Y%m%d%H')+datetime.timedelta(hours=int(frst))).strftime("%Y-%m-%d_%H:%M:%S")
intp = (datetime.datetime.strptime(date,'%Y%m%d%H')+datetime.timedelta(hours=int(step))).strftime("%Y-%m-%d_%H:%M:%S")
fstart = path+'/wrfout_d01_'+start+'*'
fintp  = path+'/wrfout_d01_'+intp+'*'
fend   = path+'/wrfout_d01_'+end+'*'
file = path+'/*'
filestart = glob.glob(fstart)
fileintp  = glob.glob(fintp)
fileend   = glob.glob(fend)
filelist  = glob.glob(file)
filelist.sort()   
rstart = np.array(np.where(np.array(filelist)==filestart))[0][0]
rintp = np.array(np.where(np.array(filelist)==fileintp))[0][0]
rend   = np.array(np.where(np.array(filelist)==fileend))[0][0]
fn = filelist[rstart:rend:rintp]
outroot = 'Users/'    
for i in fn:
    uhdr = {"header":{"discipline":0,"disciplineName":"Meteorological products","gribEdition":2,"gribLength":131858,"center":0,"centerName":"WRF OUTPUT","subcenter":0,"refTime":"2014-01-31T00:00:00.000Z","significanceOfRT":1,"significanceOfRTName":"Start of forecast","productStatus":0,"productStatusName":"Operational products","productType":1,"productTypeName":"Forecast products","productDefinitionTemplate":0,"productDefinitionTemplateName":"Analysis/forecast at horizontal level/layer at a point in time","parameterCategory":2,"parameterCategoryName":"Momentum","parameterNumber":2,"parameterNumberName":"U-component_of_wind","parameterUnit":"m.s-1","genProcessType":2,"genProcessTypeName":"Forecast","forecastTime":3,"surface1Type":103,"surface1TypeName":"Specified height level above ground","surface1Value":10,"surface2Type":255,"surface2TypeName":"Missing","surface2Value":0,"gridDefinitionTemplate":0,"gridDefinitionTemplateName":"Latitude_Longitude","numberPoints":65160,"shape":6,"shapeName":"Earth spherical with radius of 6,371,229.0 m","gridUnits":"degrees","resolution":48,"winds":"true","scanMode":0,"nx":360,"ny":181,"basicAngle":0,"subDivisions":0,"lo1":0,"la1":90,"lo2":359,"la2":-90,"dx":1,"dy":1}}

    vhdr = {"header":{"discipline":0,"disciplineName":"Meteorological products","gribEdition":2,"gribLength":131858,"center":0,"centerName":"WRF OUTPUT","subcenter":0,"refTime":"2014-01-31T00:00:00.000Z","significanceOfRT":1,"significanceOfRTName":"Start of forecast","productStatus":0,"productStatusName":"Operational products","productType":1,"productTypeName":"Forecast products","productDefinitionTemplate":0,"productDefinitionTemplateName":"Analysis/forecast at horizontal level/layer at a point in time","parameterCategory":2,"parameterCategoryName":"Momentum","parameterNumber":3,"parameterNumberName":"V-component_of_wind","parameterUnit":"m.s-1","genProcessType":2,"genProcessTypeName":"Forecast","forecastTime":3,"surface1Type":103,"surface1TypeName":"Specified height level above ground","surface1Value":10,"surface2Type":255,"surface2TypeName":"Missing","surface2Value":0,"gridDefinitionTemplate":0,"gridDefinitionTemplateName":"Latitude_Longitude","numberPoints":65160,"shape":6,"shapeName":"Earth spherical with radius of 6,371,229.0 m","gridUnits":"degrees","resolution":48,"winds":"true","scanMode":0,"nx":360,"ny":181,"basicAngle":0,"subDivisions":0,"lo1":0,"la1":90,"lo2":359,"la2":-90,"dx":1,"dy":1}}

    data = [uhdr, vhdr]
    newf = Dataset(i)
    lat = np.array(newf.variables['lat'])
    # print(fn,lat)
    lon = np.array(newf.variables['lon'])
    dys = np.diff(lat, axis = 0).mean(1)
    dy = float(dys.mean())
    dxs = np.diff(lon, axis = 1).mean(0)
    dx = float(dxs.mean())
    nx = float(lon.shape[1])
    ny = float(lat.shape[0])
    la1 = float(lat[-1, -1])
    la2 = float(lat[0, 0])
    lo1 = float(lon[0, 0])
    lo2 = float(lon[-1, -1])

    time =(newf.variables['time'])

    dates = nc.num2date(time[:],units=time.units)

    dt = pd.to_datetime(np.array(dates, dtype='datetime64[s]')).strftime("%Y%m%d%H%M%S")

    tms =pd.to_datetime(np.array(dates, dtype='datetime64[s]')).strftime("%Y-%m-%d_%H:%M:%S")
    for ti, time in enumerate(dt):

        datestr = (dt[0][:8])
        timestr = (dt[0][8:10])+'00'

        dirpath = outroot + date
        os.makedirs(dirpath, exist_ok = True)
        outpath = os.path.join(dirpath, '%s.json' % (i[-19:]))
        for u0_or_v1 in [0, 1]:

            h = data[u0_or_v1]['header']
            h['la1'] = la1
            h['la2'] = la2
            h['lo1'] = lo1
            h['lo2'] = lo2
            h['nx'] = nx
            h['ny'] = ny
            h['dx'] = dx
            h['dy'] = dy
            h['forecastTime'] = 0
            h['refTime'] = tms[0] + '.000Z'

            h['gribLength'] = 1538 + nx * ny * 2
            if u0_or_v1 == 0:
                data[u0_or_v1]['data'] = np.array(newf.variables['u10']).ravel().tolist()
            elif u0_or_v1 == 1:
                data[u0_or_v1]['data'] = np.array(newf.variables['v10']).ravel().tolist()
        if ti == 0:
            outf = open(outpath, 'w')
            json.dump(data, outf)
            outf.close()
        outf = open(outpath, 'w')
        json.dump(data, outf)
        outf.close()

上述指令碼為Linux系統下執行，執行方式如下：

python xx.py 起報時間時常間隔

舉個例子：

我的wrfout資料名稱如下：

python  convert_to_json.py 2022071000 12 06

根據你需要的模式起始時間，起報的時長（小時）以及預報的時間間隔（小時）進行自動化轉換。

python 格式轉換指令碼2

當然，這裡也準備了一個windows下的簡易指令碼，轉換出的資訊也比較簡單，

# -*- coding: utf-8 -*-
"""
Created on %(date)s

@author: %(jixianpu)s

Email : 211311040008@hhu.edu.cn

introduction : keep learning althongh walk slowly
"""
from __future__ import print_function, unicode_literals
import pandas as pd
import os
import json
import netCDF4 as nc
import numpy as np
import  datetime
from netCDF4 import Dataset
import argparse
from argparse import RawDescriptionHelpFormatter
import xarray as xr
# parser = argparse.ArgumentParser(description = """
# """, formatter_class = RawDescriptionHelpFormatter)

args = r'J:/wrf自動化/wrfout_d01_2022-07-10_01_00_00.nc'

outroot = r'D:/'

uhdr = {"header":{
                  "nx":360,
                  "ny":181,
                  "max":11,
                  }}

data = [uhdr]
newf = Dataset(args)
lat = np.array(newf.variables['lat'])
lon = np.array(newf.variables['lon'])
u10 = np.array(newf.variables['u10'])
v10 = np.array(newf.variables['v10'])

# indx = u10>1000

# u10[indx] = np.nan
# v10[indx] = np.nan

w10 = np.nanmax(np.sqrt(u10*u10+v10*v10))

dys = np.diff(lat, axis = 0).mean(1)
dy =    float(dys.mean())
print('Latitude Error:', np.abs((dy / dys) - 1).max())
print('Latitude Sum Error:', (dy / dys - 1).sum())
dxs = np.diff(lon, axis = 1).mean(0)
dx =    float(dxs.mean())
print('Longitude Error:', np.abs(dx / dxs - 1).max())
print('Longitude Sum Error:', (dx / dxs - 1).sum())

nx =    float(lon.shape[1])
ny =    float(lat.shape[0])

la1 =    float(lat[-1, -1])
la2 =   float(lat[0, 0])
lo1 =   float(lon[0, 0])
lo2 =   float(lon[-1, -1])
time =(newf.variables['time'])
dates = nc.num2date(time[:],units=time.units)
dt = pd.to_datetime(np.array(dates, dtype='datetime64[s]')).strftime("%Y%m%d%H%M%S")

ds= {
                      "nx":360,
                      "ny":181,
                      "max":11,
                      # "lo1":0,
                      # "la1":90,
                      # "lo2":359,
                      # "la2":-90,
                      # "dx":1,
                      # "dy":1,
                      # "parameterUnit":"m.s-1",
                      'data':{}
        }

ds['max']   =    float(w10)
ds['nx']    =    (nx)
ds['ny']    =    (ny)
for ti, time in enumerate(dt):
    #2012/02/07/0100Z/wind/surface/level/orthographic=-74.01,4.38,29184
    datestr = (dt[0][:8])
    timestr = (dt[0][8:10])+'00'
    print('Add "#' + datestr + '/' + timestr + 'Z/wind/surface/level/orthographic" to url to see this time')
    dirpath = os.path.join('D:', *datestr.split('/'))
    os.makedirs(dirpath, exist_ok = True)
    outpath = os.path.join(dirpath, '%s-wind-surface-level-gfs-1.0.json' % (timestr,))
    udata=u10.ravel()
    data[0]['data']=[]
    for i in range(len(udata)):

        data[0]['data'].append([
        u10.ravel().tolist()[i],
        v10.ravel().tolist()[i]])

    ds['data'] = data[0]['data']

outf = open(outpath, 'w')
json.dump(ds,outf)
outf.close()

這個指令碼正常放在編輯器裡面執行即可。

執行完結束，會在你的輸出路徑下生成一個資料夾：