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