NCL站點資料畫中國區域氣溫散點圖及分佈圖
阿新 • • 發佈:2018-12-10
參考了網上程式碼:http://bbs.06climate.com/forum.php?mod=viewthread&tid=11417 最NCL做站點資料不妨看看這個。程式碼如下:
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" begin file_path="/mnt/c/users/hong/desktop/1.txt" data=asciiread(file_path, (/4216,13/), "float") station=data(:,0) ;讀入站點號 lat=data(:,1) ;讀入緯度 lon=data(:,2) ;讀入經度 day=data(:,6) ;由於我選的就是1951年1月的,所以只有日期存在區別。 temp=data(:,7) ;NCL是從0開始計數,因此第8列索引是7。 temp@_FillValue=32766 ;在資料說明裡說了氣溫的缺測值是32766 ;由於在資料裡是整數性,需要對其進行一下轉換,其中經緯度要從60進位制轉為100進位制 lat_a=round(lat*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100 lon_a=round(lon*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100 temp_a=temp*0.1 temp_av=new(136,"float") lon_av=new(136,"float") lat_av=new(136,"float") do i=0,135 temp_av(i)=dim_avg_n(temp_a(31*i:31*i+30),0) lon_av(i)=lon_a(31*i) lat_av(i)=lat_a(31*i) end do temp_max=max(temp_av) temp_min=min(temp_av) olon=new(64,"float") ;中國經度範圍73°33′E-135°05′E,這裡我設定經度70°-140° olat=new(41,"float") ;中國緯度範圍3°51′N-53°33′N,這裡我設定緯度0°-55° data1=new((/41,64/),"float") do i=0,63 olon(i)=72.0+i end do do l=0,40 olat(l)=l+15 end do ;要設定經緯度的單位和名稱,否則會報錯 lon_av!0 = "lon" lon_av@long_name = "lon" lon_av@units = "degrees-east" lon_av&lon = lon_av lat_av!0 = "lat" lat_av@long_name = "lat" lat_av@units = "degrees_north" lat_av&lat = lat_av olon!0 = "lon" olon@long_name = "lon" olon@units = "degrees-east" olon&lon = olon olat!0 = "lat" olat@long_name = "lat" olat@units = "degrees_north" olat&lat = olat ;進行插值操作 rscan=(/20,10,5/);網上是10 5 3 我這裡用這個是為了讓每個位置有數字,因為我沒設定缺測值 R=obj_anal_ic_deprecated(lon_av,lat_av,temp_av,olon,olat,rscan, False) ;下面在開工作空間之前,設定了colormap,我直接複製過來了。 cmap = (/ \ (/ 255./255, 255./255, 255./255 /), \ ; 0 - White background. (/ 0./255 , 0./255 , 0./255 /), \ ; 1 - Black foreground. (/ 255./255, 0./255 , 0./255 /), \ ; 2 - Red. (/ 0./255 , 0./255 , 255./255 /), \ ; 3 - Blue. (/ 164./255, 244./255, 131./255 /), \ ; 4 - Ocean Blue. (/ 0./255 , 0./255 , 255./255 /), \ ; 5 - Bar 1 (/ 0./255 , 153./255, 255./255 /), \ ; 6 - Bar 2 (/ 0./255, 153./255, 153./255 /), \ ; 7 - Bar 3 (/ 0./255 , 255./255, 0./255 /), \ ; 8 - Bar 4 (/ 255./255, 255./255 , 102./255 /), \ ; 9 - Bar 5 (/ 255./255, 153./255 , 102./255 /), \ ; 10 - Bar 6 (/ 255./255, 0./255 , 255./255 /) \ ; 11 - Bar 7 /) wks=gsn_open_wks("png","example_195101") gsn_define_colormap(wks,cmap) ; res=True res@gsnAddCyclic=False ;如果設定為真,則迴圈點被加入資料,如果資料不是迴圈的,就設定為假就可以。 res@mpDataBaseVersion = "Ncarg4_1";網上的那個程式碼裡沒有這句,害我折騰了好久才明白 res@mpDataSetName="Earth..4" ;中國地圖包含在這個叫Earth..4的地相簿裡 res@mpOutlineOn=True res@mpOutlineSpecifiers=(/"China:states","Taiwan"/) res@mpGeophysicalLineThicknessF=2.0 ;這兩行是為了加粗邊界和國界線 res@mpNationalLineThicknessF=2.0 res@mpMinLatF=15.0 res@mpMaxLatF=55.0 res@mpMinLonF=72 res@mpMaxLonF=135.0 res@mpDataBaseVersion = "Ncarg4_1" res@mpAreaMaskingOn = True ;使能填充覆蓋 res@mpMaskAreaSpecifiers = (/"China:states","Taiwan"/) res@mpOceanFillColor=0 res@mpInlandWaterFillColor=0 res@cnFillOn=True;畫填充圖 res@cnLinesOn=False;不畫等值線 res@cnLineLabelsOn=False;不要等值線上的標籤 res@cnFillDrawOrder="PreDraw";先畫填充 map=gsn_csm_contour_map(wks,R,res) end
畫出的是等值線填充圖,圖如下:
散點分佈圖程式碼:
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" begin file_path="/mnt/c/users/hong/desktop/1.txt" data=asciiread(file_path, (/4216,13/), "float") station=data(:,0) ;讀入站點號 lat=data(:,1) ;讀入緯度 lon=data(:,2) ;讀入經度 day=data(:,6) ;由於我選的就是1951年1月的,所以只有日期存在區別。 temp=data(:,7) ;NCL是從0開始計數,因此第8列索引是7。 temp@_FillValue=32766 ;在資料說明裡說了氣溫的缺測值是32766 ;由於在資料裡是整數性,需要對其進行一下轉換,其中經緯度要從60進位制轉為100進位制 lat_a=round(lat*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100 lon_a=round(lon*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100 temp_a=temp*0.1 temp_av=new(136,"float") lon_av=new(136,"float") lat_av=new(136,"float") do i=0,135 temp_av(i)=dim_avg_n(temp_a(31*i:31*i+30),0) lon_av(i)=lon_a(31*i) lat_av(i)=lat_a(31*i) end do R=temp_av temp_max=max(temp_av) temp_min=min(temp_av) itv=(temp_max-temp_min)/5 arr=(/temp_min+itv,temp_min+2*itv,temp_min+3*itv,temp_min+4*itv/);我把所有溫度均勻地用4個節點分配為5份 colors = (/10,38,56,66,94/);5個水平需要5個顏色來代表 num_markers=dimsizes(arr)+1; lat_new = new((/num_markers,dimsizes(R)/),float,-999); lon_new = new((/num_markers,dimsizes(R)/),float,-999); labels = new(dimsizes(arr)+1,string);最後出現在圖下方的標籤 do i =0,num_markers-1 if(i.eq.0);第一個水平即低於0的水平 indexes=ind(R.lt.arr(0)) labels(i)="R<"+arr(0) end if if(i.eq.(num_markers-1))then;最大的一個水平即為大於26的 indexes=ind(R.ge.max(arr)) labels(i)="R>="+max(arr) end if if(i.gt.0.and.i.lt.(num_markers-1))then;中間的水平 indexes=ind(R.ge.arr(i-1).and.R.lt.arr(i)) labels(i)=arr(i-1)+"<=R<"+arr(i) end if if(.not.any(ismissing(indexes)))then;如果index裡有數,而不全是-999,那麼把lat、lon_new的前N個(在這一水平裡有N個點)換成這N個點的經緯度。 npts_range=dimsizes(indexes) lat_new(i,0:(npts_range-1))=lat_av(indexes) lon_new(i,0:(npts_range-1))=lon_av(indexes) end if delete(indexes) end do wks=gsn_open_wks("png", "scatter_example") gsn_define_colormap(wks, "WhViBlGrYeOrRe") mpres=True mpres@gsnFrame=False mpres@pmTickMarkDisplayMode="Always" mpres@mpMinLatF=15.0 mpres@mpMaxLatF=55.0 mpres@mpMinLonF=72 mpres@mpMaxLonF=135.0 mpres@tiMainString="1951 January China average temperature" mpres@mpDataBaseVersion = "Ncarg4_1";這一步和下一步必須要,否則載入中國地圖的時候會出錯(找不到地相簿) mpres@mpDataSetName="Earth..4";這步加上步再加下面那個China:state和Taiwan就可以畫出中國輪廓邊界了 mpres@mpOutlineOn=True mpres@mpOutlineSpecifiers=(/"China:states","Taiwan"/);在這個地相簿裡我們繪製中國和臺灣的區域邊界,China:state里居然沒有臺灣!不能忍,要包括進來! mpres@mpGeophysicalLineThicknessF=2.0 ;這兩行是為了加粗邊界和國界線 mpres@mpNationalLineThicknessF=2.0 ;我看網上那個最多人轉載的那個站點畫氣溫分佈等值線圖,它用的蘭伯特投影,我畫了一下,左右緯度不對稱,很奇怪,就放棄使用了。 map=gsn_csm_map(wks,mpres) gsres=True gsres@gsMarkerIndex=16 do i=0,num_markers-1 if (.not.ismissing(lat_new(i,0))) then gsres@gsMarkerColor=colors(i) gsres@gsMarkerThicknessF=i+1 gsn_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres) end if end do frame(wks); end
畫的是散點分佈圖,圖如下: