# -*- coding: utf-8 -*-
"""
Created on Mon Apr 3 09:28:07 2023
@author: PC
"""
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 11 16:54:30 2022
@author: PC
"""
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
import xarray as xr
import datetime
import os
import datetime
import metpy
from scipy.ndimage import gaussian_filter
from metpy.interpolate import cross_section
from matplotlib.pyplot import MultipleLocator
def getFullFileNameByFacotorDate(fathPath,factor,date):
# filename='''era5.{}.{}.nc'''.format(factor,date.strftime('%Y%m%d'))
# f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/catalog.xml'
filename= f'era5.{factor}.{date:%Y%m%d}.nc'
fullPath=os.path.join(fathPath,factor)
if (factor =='2_metre_dewpoint_temperature'):
# filename='''era5.{}.{}.nc'''.format('2m_dewpoint_temperature',date.strftime('%Y%m%d'))
filename= f'era5.2m_dewpoint_temperature.{date:%Y%m%d}.nc'
if (factor =='2m_Temperature'):
# filename='''era5.{}.{}.nc'''.format('2mt',date.strftime('%Y%m%d'))
filename= f'era5.2mt.{date:%Y%m%d}.nc'
return os.path.join(fullPath,filename)
# def getDate():
# f_path=r'I:\weatherExample\202108\EAR5'
# date1=datetime.datetime(2021,7,28)
# uwnd=xr.open_dataset(getFullFileNameByFacotorDate(f_path,'u_component_of_wind',date1))
# vwnd=xr.open_dataset(getFullFileNameByFacotorDate(f_path,'v_component_of_wind',date1))
# hght=xr.open_dataset(getFullFileNameByFacotorDate(f_path,'geopotential',date1))
# rh=xr.open_dataset(getFullFileNameByFacotorDate(f_path,'relative_humidity',date1))
# hght = hght/9.8
# lel_sel = 700
def get_qflux_DIV_by_ear5(uwnd,vwnd, spe_humidity):
qv_u = uwnd.u*spe_humidity/(metpy.constants.g) # g的单位为m/s**2,换算为N/kg,再换算为Pa·m2/kg,hpa.m.cm/kg,最终单层水汽通量的单位是kg/m•hPa•s
#(m.s-1. kg/hpa.m.cm)=kg/cm.hpa.s= 10**2kg/m.hpa.s
qv_v = vwnd.v*spe_humidity/(metpy.constants.g)
qv_u=qv_u*100
qv_v=qv_v*100 #最终单层水汽通量的单位是kg/m•hPa•s 如果要转为常的用10 * g/(cm.hpa.s)
a = np.sqrt(qv_u * qv_u + qv_v*qv_v)
lev= uwnd.level # 计算q*v/g,单位是kg/m•hPa•s
total_q_u = np.trapz(qv_u.q,lev,axis=1) #将单位kg/(m*s)
total_q_v = np.trapz(qv_v.q,lev,axis=1)
total_a = np.sqrt(total_q_u * total_q_u + total_q_v * total_q_v) #计算整层水汽通量 单位kg/(m*s)
#水通量散度的单位 为 g/(cm*m*s*hPa)
lev= uwnd.level
time=spe_humidity.time
lat = uwnd.latitude.data
lon = uwnd.longitude.data
lons, lats = np.meshgrid(lon, lat)
# dx 241,280 dy 240,281
div_qv = np.zeros((spe_humidity.time.shape[0],lev.shape[0],lat.shape[0],lon.shape[0]))
#(24,19,241,281)
for j in range(spe_humidity.time.shape[0]):
print(j)
for i in range(lev.shape[0]):
div_qv[j,i] = mpcalc.divergence(u = qv_u.q[j,i],v = qv_v.q[j,i])
# 单层的水汽通量散度 # 单位是kg/m2•hPa•s
# # 计算整层水汽通量散度
total_div_qv = np.trapz(div_qv[::-1],lev[::-1],axis=1)*10**5 #单位为10-5kg/(m**2*s)
#24,241,281)
div_qv_nc = xr.Dataset(
{
"mfd":(("time",'level',"latitude","longitude"), div_qv)
},
coords={
'level':lev,
"time":time,
"latitude":lat,
"longitude":lon,
}
)
div_qv_nc.attrs["long_name"] = "div-qv"
div_qv_nc.attrs["unit"] ='kg/(m**2•hPa•s)'
div_qv_nc.to_netcdf("qiv-qv.nc")
total_qv_nc = xr.Dataset(
{
"total_mfd":(("time","latitude","longitude"), total_div_qv)
},
coords={
"time":time,
"latitude":lat,
"longitude":lon,
}
)
total_qv_nc.attrs["long_name"] = "total-div-qv"
total_qv_nc.attrs["unit"] ='10-5kg/(m**2.s)'
total_qv_nc.to_netcdf("total-qiv-qv.nc")
return qv_u,qv_v,a,tota