import os
import numpy as np
import matplotlib.image as image
import matplotlib as mpl
import matplotlib.pyplot as pl
import warnings
warnings.filterwarnings('ignore')
try:
get_ipython().magic("matplotlib inline")
except:
pl.ion()
from scipy import interpolate
import datetime as dt
import glob
import json
from datetime import datetime
import pandas as pd
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, Hodograph, SkewT
from metpy.units import units
from scipy.signal import medfilt
import xarray as xr
from rpy2.robjects.packages import importr
from rpy2.robjects import r,pandas2ri
import rpy2.robjects as robjects
import rpy2.robjects.numpy2ri
rpy2.robjects.numpy2ri.activate()
pandas2ri.activate()
importr('thunder')
# import sharppy
# import sharppy.sharptab.profile as profile
# import sharppy.sharptab.interp as interp
# import sharppy.sharptab.winds as winds
# import sharppy.sharptab.utils as utils
# import sharppy.sharptab.params as params
# import sharppy.sharptab.thermo as thermo
rpy2.robjects.packages.Package as a <module 'thunder'>
with open('/home/lmathias/Documents/metpy/bufr10618/sounding.json') as json_file:
rs_data = json.load(json_file)
#with open('/home/lmathias/Desktop/SOUNDING/sounding2.json') as json_file:
#rs_data2 = json.load(json_file)
meta=rs_data[1]
year=meta[12]
print(year)
month=meta[13]
print(month)
day=meta[14]
print(day)
hour=meta[15]
print(hour)
data=rs_data[3]
data=np.array(data[2])
wmoid=data[0,3]
print(wmoid)
starthour=data[0,38]
startminute=data[0,39]
lat=data[0,41]
lon=data[0,42]
print(lat)
print(lon)
#height=data[0,43]
metdata=data[0,68:21008]
#metdata=metdata[:-13]
press=metdata[::200]/100*units.hPa
press=press.astype(float)
size = press.size
print(size)
height=metdata[1::200]
height=height.astype(float)
height=height*units.meter
height = height[0:size]
temp=(metdata[4::200]-273.15)*units.degC
temp=temp.astype(float)
print(temp.size)
dew=(metdata[5::200]-273.15)*units.degC
dew=dew.astype(float)
dew = dew[0:size]
wdir=metdata[6::200]* units.degrees
wdir = wdir.astype(int)
wdir = wdir[0:size]
print(wdir.size)
size = wdir.size
print(size)
wspd=metdata[7::200]* units('meters/second')
wspd= wspd.astype(float)
wspd = wspd[0:size]
print(wspd.size)
temp = temp[0:size]
print(temp.size)
press = press[0:size]
dew = dew[0:size]
height = height[0:size]
u, v = mpcalc.wind_components(wspd, wdir)
#q = mpcalc.specific_humidity_from_dewpoint(dew,press)
#m = mpcalc.mixing_ratio_from_specific_humidity(q)
#tempv = mpcalc.virtual_temperature(temp,m)
#tw = mpcalc.wet_bulb_temperature(press,temp,dew)
pressf=medfilt(press,9)*units.hPa
print(pressf.size)
heightf=medfilt(height,9)*units.meter
print(heightf.size)
dewf=medfilt(dew,9)*units.degC
tempf=medfilt(temp,9)*units.degC
2023 7 13 6 10618 49.69246 7.3259 78 78 78 78 78 78 78 78
# compute convective parameters
ws = np.sqrt(u**2 + v**2)*1.94384
ws = list(ws.magnitude)
pressure = list(press.magnitude)
temp2 = list(temp.magnitude)
altitude = list(height.magnitude)
dpt= list(dew.magnitude)
wd = mpcalc.wind_direction(u, v, convention='from')
wd = list(wd.magnitude)
ws = robjects.IntVector(ws)
wd = robjects.IntVector(wd)
dpt = robjects.IntVector(dpt)
altitude = robjects.IntVector(altitude)
pressure = robjects.IntVector(pressure)
temp2 = robjects.IntVector(temp2)
print(altitude)
#time2 = str(i*3+3)
robjects.r['sounding_save'](pressure,altitude,temp2,dpt,wd,ws,accuracy=1, filename='sounding_10618_'
+(datetime.strftime(datetime(year,month,day,hour,0,0),"%Y%m%d%H%M"))+'.png',
title='WMO 10618 | Idar-Oberstein | 49.6925°N 7.3259°E | '
+ datetime.strftime(datetime(year,month,day,hour,0,0),"%d.%m.%Y %H:%M UTC"),
max_speed = 30)
R[write to console]: Your display device is 10 x 6 in. It is recommended to use at least 10 x 7.5 in. plotting window or consider saving the layout into file
[1] 376 580 788 971 1176 1381 1587 1824 2052 2299 2533 2785 [13] 3044 3258 3490 3723 3954 4192 4417 4649 4856 5083 5317 5564 [25] 5793 6053 6317 6566 6825 7065 7307 7526 7752 8013 8261 8486 [37] 8712 8948 9211 9420 9624 9843 10079 10283 10494 10669 10890 11107 [49] 11304 11479 11649 11819 11994 12160 12322 12486 12659 12835 13042 13277 [61] 13452 13615 13781 13942 14110 14285 14455 14629 14796 14970 15134 15299 [73] 15465 15636 15811 15985 16171 16356
<rpy2.rinterface_lib.sexp.NULLType object at 0x7f694c7f3500> [RTYPES.NILSXP]