- 积分
- 6918
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-1-1
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 lysx 于 2017-7-25 13:08 编辑
分别用等经纬度投影和兰波托投影两种方式绘图?但是后者显示的图看上去不太正常。求助各位大侠帮忙!谢谢!
下面是代码。
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 24 07:13:57 2017
@author: 123
"""
import numpy as np
import datetime
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.font_manager import FontProperties
import os
font1=FontProperties(fname='C:\\Windows\\Fonts\\SimSun.ttc')
def initialize(resolution):
mspace=(60.0, 150.0, 60., -10.) # tutle
gspace=[0,0,0,0] # list
ospace=[0,0,0,0]
f=open('D:\\EcmwfDisplay\\zht\\systemSet.txt', 'r')
lines = f.readlines()
f.close
d1=('%d'%int(lines[0].strip()), '-%02d'%int(lines[1].strip()),
'-%02d'%int(lines[2].strip()), ' %02d'%int(lines[3].strip()))
datePredefined=datetime.datetime.strptime(''.join(d1), '%Y-%m-%d %H')
initialTime=datePredefined.strftime('%y%m%d%H.')
tt = datePredefined - datetime.timedelta(hours=36)
tt=tt.strftime('%Y%m%d%H')
newPath = ''.join(('D:\\EcmwfDisplay\\result\\', tt))
if os.path.exists(newPath):
__import__('shutil').rmtree(newPath)
tt=datePredefined.strftime('%Y%m%d%H')
newPath = ''.join(('D:\\EcmwfDisplay\\result\\', tt))
mkFolder( newPath )
validTimePredefined=int(lines[4].strip())
destPath=lines[9].strip()
ospace[0]=float(lines[5].strip())
ospace[1]=float(lines[6].strip())
ospace[2]=float(lines[7].strip())
ospace[3]=float(lines[8].strip())
ospace [ 0 ] = max ( mspace[0], ospace[0] )
ospace [ 1 ] = min ( mspace[1], ospace[1] )
ospace [ 2 ] = min ( mspace[2], ospace[2] )
ospace [ 3 ] = max ( mspace[3], ospace[3] )
gspace [ 0 ] = int( ( ospace[0] - mspace[0] ) / resolution )
gspace [ 1 ] = int ( ( ospace[1] - mspace[0] ) / resolution + 1)
gspace [ 2 ] = int ( ( mspace[2] - ospace[2] ) / resolution )
gspace [ 3 ] = int ( ( mspace[2] - ospace[3] ) / resolution + 1)
return destPath, validTimePredefined, ospace, gspace, \
datePredefined, initialTime, newPath
def mkFolder( destPath ):
isExists=os.path.exists(destPath)
if not isExists:
os.mkdir(destPath)
def dataRead(fln, space ):
fileName=''.join(fln)
if os.path.exists(fileName):
status=True
datas=np.genfromtxt(fname=fileName, skip_header=6)
datas = datas[space[2]:space[3], space[0]:space[1]]
else:
datas=-99.0
status=False
return datas, status
destPath, validTimePredefined, ospace, gspace, datePredefined, \
initialTime, newPath=initialize(0.125)
lon0 = ospace[0] + ( ospace[1] - ospace[0] )/2.0
lat0 = ospace[3] + ( ospace[2] - ospace[3] )/2.0
fig = plt.figure(figsize=(16, 10))
maps = Basemap(llcrnrlat=ospace[3], urcrnrlat=ospace[2],
llcrnrlon=ospace[0], urcrnrlon=ospace[1],
lon_0=lon0, lat_0=lat0, projection='lcc', resolution ='c')
nx = gspace[1]-gspace[0]
ny = gspace[3]-gspace[2]
lons, lats=maps.makegrid(nx,ny)
lats=lats[::-1]
xx,yy=maps(lons,lats)
datas, status=dataRead((destPath, 'LCC','\\999\\', initialTime,
'003'), gspace )
print initialTime
maps.drawmeridians(np.arange(ospace[0], ospace[1]+0.00001, 10.),
labels=[0, 0, 0, 1], fontsize=15, linewidth=0.4)
maps.drawparallels(np.arange(ospace[3], ospace[2]+0.00001, 10.),
labels=[1, 0, 0, 0], fontsize=15, linewidth=0.4)
maps.readshapefile('D:\\EcmwfDisplay\\shapefile\\bou2_4p', 'bou2_4p.shp',
linewidth=1,color='gray')
cs = maps.contourf(xx, yy, datas, [1, 3, 4, 6, 7, 10 ],
colors=['#196aff', '#00ff00', '#ffff00', '#ff5500', '#f535a0'])
maps.colorbar(cs, size='3%')
fig.savefig('D:\\cyl.png', dpi=300)
程序用到的两个文本文件:
systemSet.txt
2017
07
24
20
168
70.00
150.00
60.00
10.00
G:\ecmwf_thin\
vars.txt
LCC
两张图分别为等经纬度和兰波托投影的图。附件为文中用到的数据。
|
|