用插值处理站点数据画降水图出错

songyy

load "$NCARG_ROOT/lib/ncarg/database/Chinamap/zidingyimap.ncl"
undef("mask_data")
function mask_data(data

:numeric,fname[1]:string,outmask:logical)
local f,segments,geometry,segsDims,geomDims,geom_segIndex,geom_numSegs,segs_xyzIndex,\
segs_numPnts,numFeatures,ilat,ilon,i,j,\
shp_lat,shp_lon,data_lat,data_lon,startSegment,numSegments,seg,\
startPT,endPT,mask_val
begin
        f=addfile(fname,"r")
        segments=f->segments
        geometry=f->geometry
        segsDims=dimsizes(segments)
        geomDims=dimsizes(geometry)
        geom_segIndex=f@geom_segIndex
        geom_numSegs=f@geom_numSegs
        segs_xyzIndex=f@segs_xyzIndex
        segs_numPnts=f@segs_numPnts
        numFeatures=geomDims(0)
        data_mask=new(dimsizes(data),typeof(data),data@_FillValue)
        mask_val=1
        data_mask=0
        shp_lon=f->x
        shp_lat=f->y
        data_lat=data&$data!0$
        data_lon=data&$data!0$
        nlat=dimsizes(data_lat)
        nlon=dimsizes(data_lon)
        do i=0,numFeatures - 1
                startSegment=geometry(i,geom_segIndex)
                numSegments=geometry(i,geom_numSegs)
                do seg=startSegment,startSegment+numSegments-1
                        startPT=segments(seg,segs_xyzIndex)
                        endPT=startPT+segments(seg,segs_numPnts)-1
                        iilt_beg=ind((data_lat).lt.min(shp_lat(startPT:endPT)))
                        iilt_end=ind((data_lat).gt.max(shp_lat(startPT:endPT)))
                        iiln_beg=ind((data_lon).lt.min(shp_lon(startPT:endPT)))
                        iiln_end=ind((data_lon).gt.max(shp_lon(startPT:endPT)))
                        ilt_beg=0
                        iln_beg=0
                        ilt_end=nlat-1
                        iln_end=nlon-1
                        if(.not.any(ismissing(iilt_beg)))then
                                ilt_beg=iilt_beg(dimsizes(iilt_beg)-1)
                        end if
                        if(.not.any(ismissing(iilt_end)))then
                                ilt_end=iilt_end(0)
                        end if
                        if(.not.any(ismissing(iiln_beg)))then
                                iln_beg=iiln_beg(dimsizes(iiln_beg)-1)
                        end if
                        if(.not.any(ismissing(iiln_end)))then
                                iln_end=iiln_end(0)
                        end if
                        do ilat = ilt_beg,ilt_end
                                do ilon = iln_beg,iln_end
                                        if(data_mask(ilat,ilon).ne.mask_val.and.\
                                                gc_inout(data_lat(ilat),data_lon(ilon),\
                                                shp_lat(startPT:endPT),shp_lon(startPT:endPT)))then
                                                data_mask(ilat,ilon)=mask_val
                                        end if
                                end do
                        end do
                        delete([/iilt_beg,iilt_end,iiln_beg,iiln_end/])
                end do
        end do
        area_data=where(data_mask.eq.1,data,data@_FillValue)
        copy_VarMeta(data,area_data)
        return(area_data)
end
;
begin
        f=addfile("d:/1960-2015/var1continous-MaytoSep-rain598stn_1960_2015_daily.nc","r")
        fshape="$NCARG_ROOT/lib/ncarg/database/Chinamap/China_GuoJie_Polyline.shp"
        varnames=getfilevarnames(f)
        if(.not.any(ismissing(varnames)))then
                do i = 0,dimsizes(varnames)-1
                        printFileVarSummary(f,varnames(i))
                end do
        end if
       
        daily=f->var3
        days=f->days
        stid=f->stid
        lats=stid@lat
        lons=stid@lon
       
        yyyy=days/10000
        mmdd=days-yyyy*10000
        mm=mmdd/100
       
        nyear=max(yyyy)-min(yyyy)+1
        nday=30+31+31
        idx=ind(mmdd.ge.601.and.mmdd.le.831);截取夏季
        temp=onedtond(ndtooned(daily(:,idx)),(/dimsizes(stid),nyear,nday/))
        yr_rain=dim_sum_n(temp,2);各年夏季累计降水
        avecli=dim_avg_n(yr_rain,1);气候平均夏季年降水
       
        latN=55.
        latS=15.
        lonW=70.
        lonE=140.
       
        dlt=0.5;插值格点精度
        latnum=floattoint(180/dlt)+1;纬度格点数
        lonnum=floattoint(360/dlt);经度格点数
        latitude=latGlobeF(latnum,"lat","latitude","degrees_north");全球
        lontitude=lonGlobeF(lonnum,"lon","longtitude","degrees_east")
        lat=latitude({latS:latN});截取插值范围格点
        lon=lontitude({lonW:lonE})
        nlat=dimsizes(lat)
        nlon=dimsizes(lon)
        griddata=new((/nlat,nlon/),"float")
        griddata(:,:)=natgrid(lats,lons,avecli,lat,lon);站点数据插值
        griddata!0="lat"
        griddata&lat=lat
        griddata!1="lon"
        griddata&lon=lon
        data=mask_data(griddata,fshape,False)
        ;绘图
        wks=gsn_open_wks("png","jiangshui")
        gsn_define_colormap(wks,"MPL_YLOrBr")
        res=True
        res@gsnMaximize=True
        res@gsnDraw=False
        res@gsnFrame=False
        res@mpMinLatF=0.
        res@mpMaxLatF=55.
        res@mpMinLonF=70.
        res@mpMaxLonF=140.
        res@mpFillOn=False
        res@pmTickMarkDisplayMode="always"
        res@tmXTOn=False
        res@tmYROn=False
        res@gsnAddCyclic=False
        res@gsnCenterString=""
        res@gsnLeftString=""
        res@gsnRightString=""
        res@cnFillOn=True
        res@cnLinesOn=False
        res@tmXBLabelFontHeightF=0.01
        res@tmYLLabelFontHeightF=0.01
        res@lbTitleString="mm"
        res@lbTitlePosition="Right"
        res@lbTitleFontHeightF=.015
        res@lbTitleDirection="Across"
        plot=gsn_csm_contour_map(wks,data,res)
        plots=draw_Chinamap(wks,plot,False,True,False);绘制中国地图
       
        draw(plots)
        frame(wks)
end

dh5002

楼主解决了吗