Code-NCL-1/观测站点资料插值为格点资料.ncl at master · cycle13/Code-NCL-1 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$GEODIAG_ROOT/geodiag.ncl"

latfile="lat"
lonfile="lon"
prefile="PRE"

nlines=3056
nlat=1000
nlon=1000

lat=asciiread(latfile,(/nlines/),"float")
lon=asciiread(lonfile,(/nlines/),"float")
PRE=asciiread(prefile,(/nlines/),"float")

newlat=fspan(min(lat),max(lat),nlat)
newlat!0 = "lat"
newlat@long_name = "lat"
newlat@units = "degrees_north"
newlat&lat = newlat

newlon=fspan(min(lon),max(lon),nlon)
newlon!0 = "lon"
newlon@long_name = "lon"
newlon@units = "degrees_east"
newlon&lon = newlon

rscan = (/4, 3, 1, 0.5, 0.2, 0.05, 0.01, 0.001, 0.0001, 0.00001/)

newpre=obj_anal_ic_Wrap(lon,lat,PRE,newlon,newlat,rscan,False)
newpre!0="lat"
newpre&lat=newlat
newpre!1="lon"
newpre&lon=newlon
printVarSummary(newpre)
print(min(newpre))
print(max(newpre))

res=True

res@gsnDraw=False
res@gsnFrame=False
res@gsnSpreadColors=True
res@gsnAddCyclic=False

res@mpMinLatF=min(lat)
res@mpMaxLatF=max(lat)
res@mpMinLonF=min(lon)
res@mpMaxLonF=max(lon)

res@mpDataSetName = "Earth..4"
res@mpDataBaseVersion = "MediumRes"
res@mpGridAndLimbOn = True
res@mpGridLineDashPattern = 2


res@cnFillOn=True
res@cnLinesOn=False
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels=ispan(5,40,5)

res@pmTickMarkDisplayMode = "Always"
res@lbLabelAutoStride = True
res@gsnMaximize     = True

res@isShowProvince = True
res@isShowSouthChinaSea = False
res@isAddMask = False
res@isShowRivers = True
res@riverColor = "black"
res@riverThickness = 0.5
res@boundaryColor = "black"
res@boundaryThickness = 0.5
setup_china_map(res)


wks=gsn_open_wks("pdf","mm")
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
plot=gsn_csm_contour_map(wks,newpre,res)
attach_china_map(wks, plot)
draw(plot)