[ncl-talk] GKS contour problem

LSL lslrsgis at whu.edu.cn
Tue Apr 9 01:16:00 MDT 2019


Dear NCL Community,


I am writing for a GKS contour problem. The lines relative to this error 
might be something for projection. The complete NCL script is attached.


; res at mpLimitMode = "Corners"

; res at mpLeftCornerLatF = lat2d(nlat-1,nlon-1)

; res at mpLeftCornerLonF = lon2d(0,0)

; res at mpRightCornerLatF = lat2d(0,0)

; res at mpRightCornerLonF = lon2d(nlat-1,nlon-1)


The errors are given as follows:


GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

warning:TickMarkSetValues: Zero span between tmXBDataLeftF and 
tmXBDataRightF, turning bottom tickmarks off:[errno=1104]

warning:TickMarkSetValues: Zero span between tmXTDataLeftF and 
tmXTDataRightF, turning top tickmarks off:[errno=1104]

warning:TickMarkSetValues: Zero span between tmYLDataBottomF and 
tmYLDataTopF, turning left tickmarks off:[errno=1104]

warning:TickMarkSetValues: Zero span between tmYRDataBottomF and 
tmYRDataTopF, turning right tickmarks off:[errno=1104]

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER 51 ISSUED FROM SUBROUTINE GSVP :

--RECTANGLE DEFINITION IS INVALID

GKS ERROR NUMBER -107 ISSUED FROM SUBROUTINE GSVP :

--MAXIMUM NUMBER OF ERROR MESSAGES EXCEEDED


Any indications would be appreciated. Thanks in advance.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190409/d90b9a25/attachment.html>
-------------- next part --------------
;--LOAD OBS --
 domain ="d01"
 nlat = 62
 nlon = 91
 dir_OBS="/home/lsl/Downloads/EXPR/Scripts/NETCDF_GRID/"
; fil_OBS_tmax="SURF_CLI_CHN_TEM_DAY_GRID_0.5-MAX."+ymd+"_"+domain+".nc"
; fil_OBS_tmin="SURF_CLI_CHN_TEM_DAY_GRID_0.5-MIN."+ymd+"_"+domain+".nc"
; fil_OBS_ppt ="SURF_CLI_CHN_PRE_DAY_GRID_0.5."+ymd+"_"+domain+".nc"
  iPOS_BEG = 34

;--LOAD WRF--
 dir_WRF="/home/lsl/Downloads/DATA/wrfout_2015/Processed/"
 fil_WRF="wrfout_"+domain+"_2015_tzadjust.nc"
 fwrf  = addfile(dir_WRF+fil_WRF, "r")
 print(fwrf)

;--LOAD REF
  dir_REF = "/home/lsl/Downloads/DATA/wrfout_2015/"
  fil_REF = "wrfout_"+domain+".nc"
  fref = addfile(dir_REF+fil_REF,"r")
  lat2d = fref->XLAT(0,:,:)      ; Destination grid
  lon2d = fref->XLONG(0,:,:)
  lat2d at units ="degrees_north"
  lon2d at units ="degrees_east"


  print("loading files completed.................................")
;--Plotting options

;  pltDir  = "./PLOT_TEST/"
;
;  pltName = "OBS_to_WRF"+"_"+domain
;  pltPath = pltDir + pltName
;  pltType = "png"   ; "x11"
;  wks     = gsn_open_wks(pltType,pltPath)  
;  gsn_define_colormap(wks,"precip3_16lev")
  
;----------------------------------------------------------------------
; Set resources to zoom in on lat/lon area of interest using either
; the WRF map projection defined on the WRF file, or by setting the
; lat/lon limits ourselves.
;----------------------------------------------------------------------



;--files
 flnm  = systemfunc("cd "+dir_OBS+" ; ls "+"SURF_CLI_CHN_TEM_DAY_GRID_0.5-MAX."+"*_d01.nc")
 nflnm = dimsizes(flnm)
 print(flnm)
 print("===")
 

;  do nf=0,nflnm-1
  do nf=9,10-1
     print("Now is processing TEXT FILE "+flnm(nf)+"....................")

     ymd  = str_get_cols(flnm(nf), iPOS_BEG+0, iPOS_BEG+7)  ; yyyymmdd
     yyyy = str_get_cols(flnm(nf), iPOS_BEG+0, iPOS_BEG+3)
     mm   = str_get_cols(flnm(nf), iPOS_BEG+4, iPOS_BEG+5)
     dd   = str_get_cols(flnm(nf), iPOS_BEG+6, iPOS_BEG+7)

     

     fil_OBS_tmax="SURF_CLI_CHN_TEM_DAY_GRID_0.5-MAX."+ymd+"_"+domain+".nc"
     f = addfile(dir_OBS+fil_OBS_tmax, "r")
     
     doy = day_of_year(stringtoint(yyyy),stringtoint(mm),stringtoint(dd))
     print("doy is: "+doy)
     data_wrf = fref->TSK(0,:,:)
     ss = fwrf->DLY_TMAX(doy-1,:,:)
     data_wrf = (/ss/)
     printVarSummary(data_wrf)
     data_wrf at lat2d = lat2d
     data_wrf at lon2d = lon2d
     

     data_obs = data_wrf
     printVarSummary(data_obs)
     tt = f->PRECIP(0,:,:)
     data_obs = (/tt/)
     printVarSummary(data_obs)
     data_obs at _FillValue = -9999.
     data_obs at lat2d = lat2d
     data_obs at lon2d = lon2d



;     doy = day_of_year(stringtoint(yyyy),stringtoint(mm),stringtoint(dd))
;     print("doy is: "+doy)
;     data_wrf = fref->RAINC(0,:,:)
;     data_wrf = fwrf->DLY_TMAX(doy-1,:,:)
;     printVarSummary(data_wrf)
;     data_wrf at lat2d = lat2d
;     data_wrf at lon2d = lon2d
 
;     printVarSummary(data_obs)
;     printVarSummary(data_wrf)

  USE_WRF_PROJECTION = True
  if(USE_WRF_PROJECTION) then
  
  res               = True
  ;res = wrf_map_resources(fref,res)
  res at gsnDraw       = False  ; will panel later
  res at gsnFrame      = False
  res at gsnAddCyclic  = False  ; regional grids

  res at cnFillOn      = True
;  res at cnFillMode    = "RasterFill"
  res at cnLinesOn     = False
  res at cnLineLabelsOn= False
  res at cnMissingValFillColor   = "gray"

  res at cnLevelSelectionMode = "ExplicitLevels"
;  res at cnLevels             = (/0.1,1,2.5,5,7.5,10,12.5,15,20,25/) ; "mm/day" 
  res at cnLevels    = (/ -24, -20.,-16.,-12.,-8.,-4.,-2.,2.,4.,8.,12.,16.,20.,24/)
;  res at cnFillPalette        = (/"Snow","PaleTurquoise","PaleGreen","SeaGreen3" ,"Yellow" \   ; contour colors
;                              ,"Orange","HotPink","Red","Violet", "Purple", "Brown"/)       ; one more color than contour levels
  res at cnFillColors = (/ 4,5,6,7,8,9,10,11,12,13,14,15,16,17,18/)  ; set the colors to be used

  res at lbLabelBarOn  = False        ; will add in panel
  res at pmTickMarkDisplayMode = "Always"        ; nicer tickmarks
  res at mpDataBaseVersion     = "MediumRes"     ; better and more map outlines
  res at mpFillOn              =  False
    
  ;  print(".....................................................")
  ;  print("minmax are: "+min(XLAT)+" "+max(XLAT)+" "+min(XLON)+" "+max(XLON))



    res at mpProjection  = "LambertConformal"
    res at mpLambertParallel1F = 34.27         ; two parallels
    res at mpLambertParallel2F =108.93
; ref_lat   =  34.27,
; ref_lon   = 108.93,
    res at mpLambertMeridianF  = 105.0        ; central meridian
;    res at mpLimitMode         = "LatLon"
;    res at mpMinLatF           = min(XLAT)         ; map area
;    res at mpMaxLatF           = max(XLAT)         ; latitudes
;    res at mpMinLonF           = min(XLON)         ; and
;    res at mpMaxLonF           = max(XLON)         ; longitudes
;    print("minmax are: "+min(XLAT)+max(XLAT)+min(XLON)+max(XLON))

;	lat2d(0,0) is: 25.5428
;	lon2d(0,0) is: 95.0417
;	lat2d(nlat-1,nlon-1) is: 41.0317
;	lon2d(nlat-1,nlon-1) is: 125.792
 
;  res at mpLimitMode            = "Corners"          ; choose region of map
;  res at mpLeftCornerLatF       = XLAT(0,0)
;  res at mpLeftCornerLonF       = XLON(0,0)
;  res at mpRightCornerLatF      = XLAT(nlat-1,nlon-1)
;  res at mpRightCornerLonF      = XLON(nlat-1,nlon-1)

;  res at mpLimitMode            = "Corners"          ; choose region of map
;  res at mpLeftCornerLatF       = lat2d(nlat-1,nlon-1)
;  res at mpLeftCornerLonF       = lon2d(0,0)
;  res at mpRightCornerLatF      = lat2d(0,0)
;  res at mpRightCornerLonF      = lon2d(nlat-1,nlon-1)

   res at mpLimitMode = "Corners" 
   res at mpLeftCornerLatF = 16 
   res at mpLeftCornerLonF = 135 
   res at mpRightCornerLatF = 54 
   res at mpRightCornerLonF = 79
  
  print("lat2d(0,0) is: "+lat2d(0,0))
  print("lon2d(0,0) is: "+lon2d(0,0))
  print("lat2d(nlat-1,nlon-1) is: "+lat2d(nlat-1,nlon-1))
  print("lon2d(nlat-1,nlon-1) is: "+lon2d(nlat-1,nlon-1))

;---Change some of the defaults set by wrf_map_resources.
    res at mpGeophysicalLineColor      = "Black"
    res at mpGridLineColor             = "Black"
    res at mpLimbLineColor             = "Black"
    res at mpNationalLineColor         = "Black"
    res at mpPerimLineColor            = "Black"
    res at mpUSStateLineColor          = "Black"
    res at mpNationalLineThicknessF    = 1.0
    res at mpUSStateLineThicknessF    = 1.0
    res at mpGeophysicalLineThicknessF    = 1.0
    res at mpGridLineThicknessF    = 1.0    

    res at mpOutlineBoundarySets       = "AllBoundaries"  ; more outlines
    res at mpGeophysicalLineColor      = "black"    ; gray
    res at mpNationalLineColor         = "Black"    ; gray
    res at mpNationalLineThicknessF    = 2.0        ; 0.5
    res at mpGeophysicalLineThicknessF = 2.0        ; 0.5
    res at mpDataBaseVersion="MediumRes"
    res at mpDataSetName="Earth..4"
    res at mpOutlineSpecifiers         = (/"China:states"/)
    res at mpGridLineColor             = "Black"
    res at mpLimbLineColor             = "Black"
    res at mpPerimLineColor            = "Black"
    res at mpGridLineThicknessF    = 1.0
  else
    extra      = 0.5
    res at mpProjection  = "CylindricalEquidistant"
    res at mpMinLatF     = min(XLAT) - extra 
    res at mpMaxLatF     = max(XLAT) + extra
    res at mpMinLonF     = min(XLON) - extra
    res at mpMaxLonF     = max(XLON) + extra
    res at tiMainOffsetYF= -0.03           ; Move the title down
  end if

;--Panel Resources
  pres                    = True
  pres at gsnMaximize        = True
  pres at gsnPanelLabelBar   = True
  pres at pmLabelBarWidthF   = 0.8    ; make labelbar longer
  pres at gsnPanelMainFont   = "helvetica-bold"
  pres at lbLabelFontHeightF = 0.01

  pltres = True
  pltres at PanelPlot = True      ; Indicate these plots are to be paneled.
;--PLOT
     pltDir  = "./PLOT_TEST/"
     pltName = "OBS_to_WRF_"+ymd+"_"+domain
     pltPath = pltDir + pltName
     pltType = "png"   ; "x11"
     wks     = gsn_open_wks(pltType,pltPath)
     gsn_define_colormap(wks,"precip3_16lev")

     print(res)

     res at tiMainString = "WRF SIMU"
;      plot_wrf = gsn_csm_contour(wks,data_wrf,res)
;     plot_wrf = gsn_contour(wks,data_wrf,res)
     plot_wrf = gsn_csm_contour_map(wks,data_wrf,res)

     res at tiMainString = "OBS GRID"
;      plot_obs = gsn_csm_contour(wks,data_obs,res)
;     plot_obs = gsn_contour(wks,data_obs,res)
     plot_obs = gsn_csm_contour_map(wks,data_obs,res)

     if(USE_WRF_PROJECTION) then
        pres at gsnPanelMainString = ymd+"OBS vs WRF"+" "+domain + " WRF map projection"
     else
        pres at gsnPanelMainString = ymd+"OBS vs WRF"+" "+domain + " Cylindrical equidistant projection"
     end if

;--Draw both plots in a single panel

     gsn_panel(wks,(/plot_wrf,plot_obs/),(/1,2/),pres)

  end do


More information about the ncl-talk mailing list