trouble regridding data

From: Oleksandra Hararuk <shurahararuk_at_nyahnyahspammersnyahnyah>
Date: Wed May 23 2012 - 14:51:46 MDT

Hello,

I am trying to regrid data from high resolution to lower resolution,
however I keep getting an error "Error in "cremapbin": could not map global
lat array into grid array". There was a thread about this in September,
however, my data are global and I get same error when I use
area_hi2lores_Wrap.
Please help.
My code is minimally modified example from the ncl webpage (regrid_6.ncl)
and I attach the input file as well:

;*************************************************************
; regrid_6.ncl
;
; Concepts illustrated:
; - Reading data of type "short" and converting it to float
; - Reordering North->South data so it is South->North
; - Interpolating from a fixed grid to a lower resolution grid using
conservative remapping.
; - Computing global areal mean values by generating the weights
; - Creating a color map using named colors
; - Drawing color-filled contours over a cylindrical equidistant map
; - Paneling three plots vertically on a page
; - Creating a netCDF file
; - Using "systemfunc" to execute a UNIX command
; - Using "system" to execute a UNIX command
; - Using "systemfunc" to get the current date
;*************************************************************

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
;*************************************************************
; User specified options
;*************************************************************
  PLOT = True
  netCDF = True

  pltDir = "./" ; only used if PLOT=True, otherwise ignored
  pltName = "regrid"
  pltType = "pdf" ; x11, ps, pdf, eps, png [v5.2.0]

  ncDir = "./" ; only used if netCDF=True, otherwise ignored
  ncFil = "GPCP_2x3"
;*************************************************************
; read and reorder the data
;*************************************************************
   diri = "./"
   fili = "7yearaveNPP.nc" ; (lat, lon) => (180,360)

   f = addfile(diri+fili, "r")
   x = f->NPP; x = short2flt(f->NPP) ;
(time,lat,lon), original order N->S
   x = x(::-1,:) ; Make s->N as required by
"area_conserve_remap"

 ; printVarSummary(x)
 ; printMinMax(x, True)
;*************************************************************
; perform conservative remapping to two different grid resolutions
;*************************************************************
   opt = False
      NLAT2x3= 64;90 ; RES = "2x3"
   MLON2x3= 128;120
   ;LAT2x3 = latGlobeFo(NLAT2x3, "LAT", "latitude" , "degrees_north") ;
89S -> 89N
  ; LON2x3 = lonGlobeFo(MLON2x3, "LON", "longitude", "degrees_east" ) ;
1.5E->358.5E
  printVarSummary(LAT2x3)
printVarSummary(LON2x3)
print(LAT2x3)
   X2x3 = area_conserve_remap_Wrap (x&lon, x&lat, x ,LON2x3, LAT2x3, opt)
printVarSummary(X2x3)

; ---
   NLAT5x5= 36 ; RES = "5x5"
   MLON5x5= 72
   LAT5x5 = latGlobeFo(NLAT5x5, "LAT", "latitude" , "degrees_north") ;
87.5S -> 87.5N
   LON5x5 = lonGlobeFo(MLON5x5, "LON", "longitude", "degrees_east" ) ;
2.5E->357.5E
   X5x5 = area_conserve_remap_Wrap (x&lon, x&lat, x ,LON5x5, LAT5x5, opt)

   printVarSummary(X2x3)
   printMinMax(X2x3, True)

   printVarSummary(X5x5)
   printMinMax(X5x5, True)

;*************************************************************
; For illustration, compute the global means of input and output grids
;*************************************************************
   gwi = latRegWgt(x&lat , "double", 0) ; input lat weights
   gwo2x3 = latRegWgt(LAT2x3, "double", 0) ; output lat weights
2x3
   gwo5x5 = latRegWgt(LAT5x5, "double", 0) ; output lat weights
5x5

   xAvg1x1 = wgt_areaave (x , gwi , 1.0, 0)
   xAvg2x3 = wgt_areaave (X2x3, gwo2x3, 1.0, 0)
   xAvg5x5 = wgt_areaave (X5x5, gwo5x5, 1.0, 0)

   xAvgDiff2x3 = xAvg2x3-xAvg1x1
   xAvgDiff5x5 = xAvg5x5-xAvg1x1

   print(xAvg1x1+" "+xAvg2x3+" "+xAvg5x5+" "+xAvgDiff2x3+"
"+xAvgDiff5x5)

   eps = 0.001
   if (max(abs(xAvgDiff2x3)).lt.eps) then
       print("area_conserve_remap: 1x1 => 2x3: SUCCESS")
   else
       print("area_conserve_remap: 1x1 => 2x3: FAIL: maxDiff
="+max(abs(xAvgDiff2x3)))
   end if
   if (max(abs(xAvgDiff5x5)).lt.eps) then
       print("area_conserve_remap: 1x1 => 5x5: SUCCESS")
   else
       print("area_conserve_remap: 1x1 => 5x5: FAIL: maxDiff
="+max(abs(xAvgDiff5x5)))
   end if

;***************************************************************
; Create plot
;***************************************************************
   if (PLOT) then
       wks = gsn_open_wks(pltType, pltDir+pltName)
       colors = (/"white","black","Snow" \ ; "WhiteSmoke" \
               ,"PaleTurquoise","PaleGreen","SeaGreen3" ,"Yellow" \
               ,"Orange","HotPink","Red","Violet", "Purple", "Brown"/)
       gsn_define_colormap(wks, colors) ; generate new color
map

       plot = new ( 3 , "graphic")

       res = True ; plot mods desired
       res@gsnDraw = False ; don't draw
       res@gsnFrame = False ; don't advance frame

       res@cnFillOn = True ; turn on color fill
       res@cnLinesOn = False ; turn of contour lines
      ;res@cnFillMode = "CellFill" ; Cell Mode
       res@cnFillMode = "RasterFill" ; Raster Mode
       res@cnLinesOn = False ; Turn off contour lines
       res@cnLineLabelsOn = False ; Turn off contour lines
       res@cnLevelSelectionMode = "ExplicitLevels"
       res@cnLevels = (/0.1,1,2.5,5,10,15,20,25,50,75/) ;
"mm/day"
       res@cnMissingValFillPattern = 0 ; make 'missing' black
       res@cnMissingValFillColor = "black"

       res@lbLabelBarOn = False ; turn off individual cb's

       res@mpCenterLonF = 210.
       res@mpFillOn = False

      ;res@mpShapeMode = "FreeAspect"
      ;res@vpWidthF = 0.8
      ;res@vpHeightF = 0.4

       nt = 0

       res@gsnLeftString = "1x1"
       res@gsnCenterString = "Areal Mean="+sprintf("%6.5f", xAvg1x1(nt))
       plot(0) = gsn_csm_contour_map_ce(wks,x(nt,:,:), res)

       res@gsnLeftString = "2x3"
       res@gsnCenterString = "Areal Mean="+sprintf("%6.5f", xAvg2x3(nt))
       plot(1) = gsn_csm_contour_map_ce(wks,X2x3(nt,:,:), res)

       res@gsnLeftString = "5x5"
       res@gsnCenterString = "Areal Mean="+sprintf("%6.5f", xAvg5x5(nt))
       plot(2) = gsn_csm_contour_map_ce(wks,X5x5(nt,:,:), res)

       resP = True
       resP@gsnMaximize = True ; make ps/eps/pdf
large [no effect x11]
     ;;resP@gsnPaperOrientation = "Portrait" ; force portrait
       resP@gsnPanelLabelBar = True ; add common colorbar
       resP@lbLabelAutoStride = True
       resP@lbLabelFontHeightF = 0.0175 ; change font size
       resP@txString = "Conservative Remap: Fixed-to-Fixed"
       gsn_panel(wks,plot,(/3,1/),resP) ; now draw as one plot
   end if ; PLOT

;************************************************
; Create netCDF ? ... Only do 2x3
; Recommend to always create a 'time' dimension
; Save only the interpolated CMORPH (uncomment to sabe "COMB")
;************************************************
   if (netCDF) then
       nline = inttochar(10)

       dimx = dimsizes(X2x3)
       ntim = dimx(0)
       nlat = dimx(1)
       mlon = dimx(2)

       time = f->time

       globeAtt = 1
       globeAtt@title = "GPCP: 1x1 interpolated to a 2x3 grid"
       globeAtt@source_file = fili
       globeAtt@creation_date= systemfunc ("date" )

       NCFILE = ncDir + ncFil +".nc"
       system ("/bin/rm -f " + NCFILE) ; remove any pre-exist file

       ncdf = addfile(NCFILE,"c")
       fileattdef( ncdf, globeAtt ) ; create the global [file]
attributes

       filedimdef(ncdf,"time",-1,True) ; make time and UNLIMITED
dimension
                                           ; recommended for most
applications
       ncdf->PRC = X2x3

   end if ; netCDF

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Wed May 23 14:52:26 2012

This archive was generated by hypermail 2.1.8 : Fri May 25 2012 - 08:35:50 MDT