怎么用python画天气图_用Python绘制天气图

痛定思痛。 2022-11-01 10:51 258阅读 0赞

\#!/bin/env python\#coding=utf-8import codecsimport numpy as npimport datetimeimport reimport pandas as pdimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom metpy.calc import wind\_componentsfrom metpy.plots import add\_metpy\_logo, current\_weather, sky\_cover, StationPlot, wx\_code\_mapfrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterfrom metpy.units import unitsimport matplotlib.pyplot as pltfrom matplotlib.patches import Polygonimport cartopy.io.shapereader as shpreaderdef func\_return\_self(x):return xdef split(str\_in, sep=None, convert\_func=None):ss = str\_in.split(sep)nn = len(ss)rst = \[None\]\*nnif convert\_func is None:cf = \[func\_return\_self\] \* nnelse:cf = \[\_ for \_ in convert\_func\]\# print(nn, len(cf))for i in range(nn):xcf = func\_return\_self if cf\[i\] is None else cf\[i\]rst\[i\] = xcf(ss\[i\])return rstdef d1\_parse(filename,encoding="gb2312"):"""读micaps第1类数据文件到内存:return:"""with codecs.open(filename, mode='r', encoding=encoding) as fid:contents = fid.readlines()dataflag, style, title = split(contents\[0\],None,None)nn = len(contents)nx = int((nn-2)/2)data = \{"stid":\[None\]\*nx,"lon":\[None\]\*nx,"lat":\[None\]\*nx,"height":\[None\]\*nx,"stlevel":\[None\]\*nx,"cloud":\[None\]\*nx,"wd":\[None\]\*nx,"ws":\[None\]\*nx,"psfc":\[None\]\*nx,"pd3":\[None\]\*nx,"pw1":\[None\]\*nx,"pw2":\[None\]\*nx,"prc6":\[None\]\*nx,"lc1":\[None\]\*nx,"lc2":\[None\]\*nx,"lc3":\[None\]\*nx,"dp":\[None\]\*nx,"njd":\[None\]\*nx,"cw":\[None\]\*nx,"tc":\[None\]\*nx,"mc":\[None\]\*nx,"hc":\[None\]\*nx,"flag1":\[None\]\*nx,"flag2":\[None\]\*nx,"dtc24":\[None\]\*nx,"dp24":\[None\]\*nx\}data = pd.DataFrame(data)j=0for i in range(2,nn,2):iline1 = contents\[i\]iline2 = contents\[i+1\]iline = " ".join(\[iline1,iline2\])\#for iline in contents\[2:\]:data.stid\[j\], data.lon\[j\], data.lat\[j\], data.height\[j\], data.stlevel\[j\],data.cloud\[j\],data.wd\[j\],data.ws\[j\],data.psfc\[j\],data.pd3\[j\], data.pw1\[j\], data.pw2\[j\], \\data.prc6\[j\], data.lc1\[j\],data.lc2\[j\],data.lc3\[j\], data.dp\[j\], data.njd\[j\], data.cw\[j\], data.tc\[j\], data.mc\[j\], data.hc\[j\], data.flag1\[j\], data.flag2\[j\], data.dtc24\[j\],data.dp24\[j\] \\=split(iline,None,\[int, float,float,float,int,float,float,float,float,float,int,int,\\float,float,float,float,float,float,int, float,float, float,int, int,float,float\])j+=1\#data.append(\[stid, lon, lat, height, stlevel,cloud,wd,ws,psfc,pd3, pw1, pw2,\\\#prc6, lc1,lc2,lc3, dp, njd, cw, tc, mc, hc, flag1, flag2, dtc24,dp24\])\#print(data)return dataflag, style, title,datadef d4\_parse(filename,encoding="gb2312"):"""读micaps第4类数据文件到内存:return:"""try:file\_object = codecs.open(filename, mode='r', encoding=encoding)all\_the\_text = file\_object.read().strip()file\_object.close()contents = re.split('\[\\s\]+', all\_the\_text)if len(contents) < 23:returndataflag = contents\[0\].strip()style = contents\[1\].strip()title = contents\[2\].strip()yy = int(contents\[3\].strip())mm = int(contents\[4\].strip())dd = int(contents\[5\].strip())hh = int(contents\[6\].strip())forehh = int(contents\[7\].strip())level = contents\[8\].strip()deltalon = float(contents\[9\].strip())deltalat = float(contents\[10\].strip())beginlon = float(contents\[11\].strip())endlon = float(contents\[12\].strip())beginlat = float(contents\[13\].strip())endlat = float(contents\[14\].strip())sumlon = int(contents\[15\].strip())sumlat = int(contents\[16\].strip())distance = float(contents\[17\].strip())min = float(contents\[18\].strip())max = float(contents\[19\].strip())def1 = contents\[20\].strip()def2 = contents\[21\].strip()x = np.arange(beginlon, endlon + 0.9\*deltalon, deltalon)y = np.arange(beginlat, endlat + 0.9\*deltalat, deltalat)X, Y = np.meshgrid(x, y)if dataflag == 'diamond' and style == '4':begin = 22Z = np.zeros((sumlat, sumlon))for i in range(sumlon):for j in range(sumlat):Z\[j, i\] = float(contents\[begin + j \* sumlon + i\])return dataflag, style,title,yy,mm,dd,hh,forehh,level,X,Y,Z,distance,min,max,def1,def2except Exception as err:print(u'【\{0\}】\{1\}-\{2\}'.format(filename, err, datetime.datetime.now()))if \_\_name\_\_ == "\_\_main\_\_":dataflag, style, title, yy, mm, dd, hh, forehh, level, X, Y, p0\_data, distance, dmin, dmax, def1, def2 = d4\_parse("/home/hliang/share/share/plotest/20190126/p0/19012602.000")station\_dataflag, staiton\_style, station\_title, station\_data = d1\_parse("/home/hliang/share/share/plotest/20190126/plot/19012108.000")station\_data = station\_data.replace(9999,np.nan)station\_data = station\_data\[station\_data.stlevel==2\]\#.dropna()lat\_max=50lat\_min=15lon\_max=130lon\_min=100proj = ccrs.LambertCylindrical(central\_longitude=(lon\_max+lon\_min)/2.0)\#ccrs.LambertConformal(central\_longitude=(lon\_max+lon\_min)/2.0, central\_latitude=(lat\_max+lat\_min)/2.0,standard\_parallels=\[(lat\_max+lat\_min)/2.0\])fig = plt.figure(figsize=(10, 10))ax = fig.add\_subplot(1, 1, 1, projection=proj)ax.set\_extent(\[lon\_min,lon\_max,lat\_min,lat\_max\],proj)ax.set\_xticks(range(lon\_min,lon\_max+2,10), crs=proj)ax.set\_yticks(range(lat\_min,lat\_max+2,10), crs=proj)shp = shpreader.Reader('/home/hliang/share/share/plotest/shp\_resources/CHN\_adm/CHN\_adm1.shp')ax.add\_geometries(shp.geometries(), crs=proj, facecolor="none", edgecolor="grey", lw=0.8)contours = ax.contour(X, Y, p0\_data,levels=np.arange(dmin,dmax+0.9\*distance,distance), colors="black",transform=proj,linewidths=0.5)plt.clabel(contours, fontsize=8,fmt='%.1f')stationplot = StationPlot(ax, station\_data\['lon'\].values, station\_data\['lat'\].values, clip\_on=True,transform=proj, fontsize=6)u,v = wind\_components(station\_data\['ws'\].values \* units('m/s'),station\_data\['wd'\].values \* units.degree)stationplot.plot\_parameter('NW', station\_data\['tc'\], color='red')stationplot.plot\_barb(u, v,sizes=dict(emptybarb=0.15))cw = station\_data.cw.replace(np.nan,0).astype(np.int)pw1 = station\_data.pw1.replace(np.nan,0).astype(np.int)stationplot.plot\_symbol('W', cw, current\_weather,color='fuchsia')stationplot.plot\_symbol('SE', pw1, current\_weather,color='y')'''for record,state in zip(shp.records(), shp.geometries()):name = record.attributes\["NAME\_1"\]if name=="Jiangsu":ax.add\_geometries(state,crs=proj,facecolor="none",edgecolor="grey",lw=0.5)'''plt.savefig('./fig1.png',dpi=300)\#plt.show()