ncl读取并绘制PIOMAS逐日的海冰厚度和海冰运动矢量资料

爱睡觉的青团

begin

;=========read PIOMAS ice thickness + ice velocity ==========================
    nlat = 120
    nlon = 360
;标量场经纬度信息 ---------------      
    files = "PIOMAS/Sea_ice_thickness/grid.dat"
    latlon2d = asciiread(files,(/2,nlat,nlon/),"float")     
    lon2d    = latlon2d(0,:,:)   
    lat2d    = latlon2d(1,:,:)   
    lat2d@units = "degrees_north"
    lon2d@units = "degrees_east"
    printVarSummary(lon2d)
    printVarSummary(lat2d)
    printMinMax(lat2d, True)  
    delete(latlon2d)

;矢量场经纬度信息
    files1 = "PIOMAS/Sea_ice_velocity/grid.dat.pop"
    data = asciiread(files1,(/7,nlat,nlon/),"float")     
    lat2d_  = data(0,:,:)   
    lon2d_ = data(1,:,:)
    lat2d_@units = "degrees_north"
    lon2d_@units = "degrees_east"
    printVarSummary(lon2d_)
    printVarSummary(lat2d_)
    printMinMax(lat2d_, False)
    printMinMax(lon2d_,False)  

    ; HTN, HTE, HUS, HUW are lengths of the 4 sides of a grid cell in km
    ; c HTN, HTE are lengths of the northern and eastern sides of a scaler grid cell in km, HTN*HTE is the area of a scaler grid cell in km**2 and can be used to calculate sea ice volume and volumes of other variables
    ; c HUS, HUW are lengths of the southern and western sides of a vector grid cell in km  
    ; HTN = data(2,:,:)
    ; HTE = data(3,:,:)
    ; HUS = data(4,:,:)
    ; HUW = data(5,:,:)
    ; HTN@units  = "km"
    ; HTE@units  = "km"
    ; HUS@units  = "km"
    ; HUW@units  = "km"
    ; printMinMax(HTN,False)
    ; printMinMax(HTE,False)
    ; printMinMax(HUS,False)
    ; printMinMax(HUW,False)

    ;angle is the angle between latitude line and  grid cell x-coordinate line, needed for plotting vecto
    angle = data(6,:,:)
    angle@units= "degree"
    printVarSummary(angle)
    printMinMax(angle,False)
    angle =  angle*3.14159265/180
    delete(data)

;掩码（陆地/海洋）ocean levels > 0, land <= 0 (文件读出来有问题，直接查看txt有120*721，但是读取出来只有20000多个数据)

    ; files2 = "PIOMAS/io.dat_360_120.output"
    ; kmt = toint(asciiread(files2, (/nlat, nlon/), "integer"))
    ; kmt@units = "land/ocean_mask"
    ; kmt@info  = "ocean > 0, land <= 0"
    ; printVarSummary(kmt)
    ; printMinMax(kmt, True)

;read ice thickness + ice velocity
    fils = systemfunc("ls PIOMAS/Sea_ice_thickness/hiday.H*")
    fils1 = systemfunc("ls PIOMAS/Sea_ice_velocity/uiday.H*")

    nyear = 2020-1979+1
    nday = 365
    thickness = new((/nyear,nday,nlat,nlon/), "float",-999.)
    uv = new((/2,nyear,nday,nlat,nlon/), "float",-999.)

    do i = 0, 0;nyear-1
        thickness(i,:,:,:) = fbindirread(fils(i), 0, (/nday,nlat,nlon/), "float")  ;m   
        uv(:,i,:,:,:) = fbindirread(fils1(i), 0, (/2,nday,nlat,nlon/), "float") ;m/s     
    end do

    thickness!2 = "lat2d"
    thickness!3 = "lon2d"

    thickness@lat2d = lat2d
    thickness@lon2d = lon2d
    printVarSummary(thickness)

    uv!3 = "lat2d"
    uv!4 = "lon2d"
    uv@lat2d = lat2d_
    uv@lon2d = lon2d_
    printVarSummary(uv)
    printMinMax(uv,False)
    u = uv(0,:,:,:,:)
    v = uv(1,:,:,:,:)
    delete(uv)

;uv rot 角度旋转
    urot = u
    vrot = v
    do i = 0, 0;nyear-1
        do j = 0, 0;nday-1
            urot(i,j,:,:) = u(i,j,:,:) * cos(angle) + v(i,j,:,:) * sin(angle)
            vrot(i,j,:,:) = -u(i,j,:,:) * sin(angle) + v(i,j,:,:) * cos(angle)
        end do
    end do
    printVarSummary(urot)
    printMinMax(urot(0,0,:,:),False)
    delete(u)
    delete(v)

;plot example
    wks = gsn_open_wks("pdf","5sample")            ; send graphics to PNG file
    gsn_define_colormap(wks,"BlueDarkRed18")

    res                     = True              ; Plot modes desired.
    res@gsnMaximize         = True              ; Maximize plot
    res@gsnAddCyclic = True
    res@gsnDraw = False
    res@gsnFrame = False
    res@gsnPolar   = "NH"                          ; specify the hemisphere
    res@mpMinLatF  =  50
    res@cnFillOn            = True              ; color plot desired
    res@cnLinesOn           = False             ; turn off contour lines
    res@cnLevelSelectionMode= "ExplicitLevels"

    ; res@cnLevels = (/"8","7","6","5","4","3","2","1"/)
    ; res@cnFillColors =  (/2,3,4,5,6,7,8,9,0/)
    res@cnLevels = (/"0.1","1","2","3","4","5","6","7"/)
    res@cnFillColors =  (/0,9,8,7,6,5,4,3,2/)


    ; 绘制海冰速度矢量
    vres = True
    vres@gsnDraw = False
    vres@gsnFrame = False
    vres@vcRefMagnitudeF = 0.2          ; 参考矢量大小
    vres@vcRefLengthF = 0.02            ; 参考矢量长度
    vres@vcMinDistanceF = 0.025          ; 矢量密度

    plot = gsn_csm_contour_map_polar(wks,thickness(0,0,:,:),res)   
    plot1 = gsn_csm_vector(wks, urot(0,0,:,:), vrot(0,0,:,:), vres)
    overlay(plot, plot1)


    resP = True
    gsn_panel(wks,plot,(/1,2/),resP)
end