python轉換wrf輸出的數(shù)據(jù)為網(wǎng)頁可視化json格式
更新時間：2022年09月26日 14:22:30 作者：oceanography-Rookie
這篇文章主要介紹了python轉換wrf輸出的數(shù)據(jù)為網(wǎng)頁可視化json格式，文章圍繞主題展開詳細的內容介紹，具有一定的參考價值，需要的小伙伴可以參考一下
前言
- 一般網(wǎng)頁可視化風場中的數(shù)據(jù)都是json格式，而如果我們希望將wrf模式模擬輸出的風場數(shù)據(jù)在網(wǎng)頁中進行展示，這就需要先將wrfoutput數(shù)據(jù)轉換為網(wǎng)頁可以識別的json格式。
- 這里主要需要用到json庫，主要的實現(xiàn)方式就是將讀取的風場風量U，V轉換為字典并存到json文件中
- 同時，由于wrf模擬的數(shù)據(jù)一般是非等間距的網(wǎng)格，需要先將數(shù)據(jù)進行插值，插值到等間距的網(wǎng)格，這里可以通過NCL的函數(shù)
rcm2rgrid_Wrap實現(xiàn)
舉個例子，將模式中設置為蘭伯特投影的網(wǎng)格：
插值為等間距網(wǎng)格：
主要的編程分為兩部分：
- 第一部分通過NCL腳本將wrfout數(shù)據(jù)轉換為等間距網(wǎng)格，并導出為netcdf格式；
- 第二部分通過python腳本將第一步導出的nc格式進行轉換，并保存輸出為json格式。
NCL插值腳本1
需要修改的就是路徑和變量，我下面展示腳本不僅有風場數(shù)據(jù)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上述腳本的缺點在于只能基于模式模擬的經緯度區(qū)域進行插值，意思就是說他的經緯度區(qū)域是固定的那么大
NCL插值腳本2
NCL還有一個函數(shù)可以實現(xiàn)上述過程，就是ESMF_regrid，該函數(shù)的優(yōu)點在于可以實現(xiàn)任意經緯度范圍的插值，但是不足在于對于存在高度層的變量，暫時無法進行高度層的數(shù)據(jù)讀取。
（也可能我水平有限不知道。。。。）這里也附上腳本：
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
endPS:運行該腳本會生成四個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 數(shù)據(jù),并生成對應文件名的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/'#只能是已經存在的文件目錄且有數(shù)據(jù)才可以進行讀取
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系統(tǒng)下運行，運行方式如下：
python xx.py 起報時間 時常 間隔
舉個例子：
我的wrfout數(shù)據(jù)名稱如下：
python convert_to_json.py 2022071000 12 06
根據(jù)你需要的模式起始時間，起報的時長（小時）以及預報的時間間隔（小時）進行自動化轉換。
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()這個腳本正常放在編輯器里面運行即可。
運行完結束，會在你的輸出路徑下生成一個文件夾：
里面有個json數(shù)據(jù)：
數(shù)據(jù)信息比較簡單，只有nx（經度的大?。?，ny（緯度的大?。┮约白畲笾担?/strong>
ok，以上就是完整的過程，最終將得到的json數(shù)據(jù)通過.js腳本運行就可以部署到網(wǎng)頁上了，簡單試了一下，大概如下圖所示，可以根據(jù)需要自行更改設置：
到此這篇關于python轉換wrf輸出的數(shù)據(jù)為網(wǎng)頁可視化json格式的文章就介紹到這了,更多相關python可視化json格式內容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關文章希望大家以后多多支持腳本之家！
相關文章
這篇文章主要介紹了python數(shù)據(jù)分析之線性回歸選擇基金，文章基于前幾篇的文章內容展開詳細的基金的趨勢分析，感興趣的小伙伴可以參考一下2022-05-05
python mongo 向數(shù)據(jù)中的數(shù)組類型新增數(shù)據(jù)操作
這篇文章主要介紹了python mongo 向數(shù)據(jù)中的數(shù)組類型新增數(shù)據(jù)操作，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧2020-12-12
今天小編就為大家分享一篇python3實現(xiàn)繪制二維點圖，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧2019-12-12
這篇文章主要為大家詳細介紹了tensorflow實現(xiàn)KNN識別MNIST，具有一定的參考價值，感興趣的小伙伴們可以參考一下2018-03-03
使用python+pygame實現(xiàn)中秋節(jié)動畫效果
馬上就要中秋節(jié)了,使用python可以實現(xiàn)中秋節(jié)動畫效果,包括月亮、兔子和煙花嗎？當然是可以的,那該如何實現(xiàn)呢？這篇文章我們主要使用pygame來實現(xiàn),文中有詳細的代碼示例供大家參考,需要的朋友可以參考下2023-09-09
本文適合編程語言零基礎的初學者，有打算轉行學習python的可以添加關注，后續(xù)小編將會把自己轉行學pyhton語言以來的所有筆記，和工作中總結的一些開發(fā)經驗分享給大家2021-09-09
Python切換pip源兩種方法(解決pip?install慢)
這篇文章主要給大家介紹了關于Python切換pip源兩種方法(解決pip?install慢),我總結的這幾種更換pip源的常用方式,希望可以幫助您成功配置國內源,解決安裝Python包速度慢的問題,需要的朋友可以參考下2023-11-11
這篇文章主要介紹了Pyinstaller加密打包應用的示例代碼，代碼簡單易懂，非常不錯，對大家的學習或工作具有一定的參考借鑒價值,需要的朋友可以參考下2020-06-06
這篇文章主要介紹了http請求 request失敗自動重新嘗試代碼示例，小編覺得還是挺不錯的，具有一定借鑒價值，需要的朋友可以參考下2018-01-01
Python編程實現(xiàn)雙擊更新所有已安裝python模塊的方法
這篇文章主要介紹了Python編程實現(xiàn)雙擊更新所有已安裝python模塊的方法,涉及Python針對模塊操作命令的相關封裝與調用技巧,需要的朋友可以參考下2017-06-06