[ncl-talk] Reading PRISM 4km data onto regular lat/lon grid

[Scott Capps]

Greetings,

I am trying to read in PRISM 4km data in bil format and output on the 
regular lat/lon grid, saving as a NetCDF file.  I eventually want to 
mask the data using a shapefile and want to keep it on a regular lat/lon 
grid.

The script works fine but I am concerned about the issues below where 
the dx and dy do not equate to 0.04166667 degrees (See lines enclosed 
by: *************).  I can provide the PRISM data when needed.  Any help 
is appreciated.

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
; USAGE
; ncl 'dir_in="./"' 'fil_in="PRISM_ppt_stable_4kmD2_19820210_bil.bil"' 
'dir_out="./"' 'file_out="prism_socalif_ppt_4km_02101982"' 
prism_to_netcdf.ncl
;
begin
;----------------------------------------------------------------------------
; FROM THE HEADER FILE:
;
; BYTEORDER      I
; LAYOUT         BIL
; NROWS          621
; NCOLS          1405
; NBANDS         1
; NBITS          32
; BANDROWBYTES   5620
; TOTALROWBYTES  56201
; PIXELTYPE      FLOAT
; ULXMAP         -125
; ULYMAP         49.9166666666687
; XDIM           0.04166666666667
; YDIM           0.04166666666667
; NODATA         -9999
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
   NBANDS = 1
   NROWS  = 621
   NCOLS  = 1405
   ORDER  = "I"
   TYPE   = "float"
   NODATA = -9999.0

   ULXMAP = -125.0d0
   ULYMAP = 49.9166666666687d0

   XDIM   = 0.04166666666667d0
   YDIM   = 0.04166666666667d0

   diri  = dir_in
   fili  = fil_in

   if (ORDER.ne."I") then
       setfileoption("bin","ReadByteOrder","LittleEndian") ; perform 
byte swapping
   end if

   pltType = "ps"
   pltDir  = "./"
   ;
   ; I am on a littleendian machine
   ;if (isbigendian() ) then
   ;  print("bigendian")
   ;else
   ;  print("littleendian")
   ;end if
   ;
   do NB=0,NBANDS-1
      dat   = fbindirread(diri+fili,NB, (/NROWS,NCOLS/), TYPE)
      dat at _FillValue = NODATA
      dat at long_name = "Daily Total Precipitation"
      dat at units     = "mm"
      dat = dat(::-1,:)
   end do
   ; 4km resolution (0.04166667x0.04166667 decimal degrees)
   lat=fspan(24.0625d0,49.9375d0,NROWS)  ; Bounded by 24.0625 to 49.9375 
Degrees North
   lon=fspan(-125.0208333d0,-66.4791667d0,NCOLS) ; Bounded by 
-125.0208333 to -66.4791667 Degrees East
   ; ***********************************
   ; Should be intervals of 0.04166667 (or, close to that at least)
   ;   instead it has unequal intervals on the order of: 0.04173387096774661
   ; ***********************************
   print(lat(1:NROWS-1)-lat(0:NROWS-2))
   ; ***********************************
   ; Should be intervals of 0.04166667
   ;   instead it has unequal intervals on the order of: 0.0416963437
   ; ***********************************
   print(lon(1:NCOLS-1)-lon(0:NCOLS-2))
   ;
   lat!0="lat"
   lat at units="degrees_north"
   lat&lat=lat
   lon!0="lon"
   lon at units="degrees_east"
   lon&lon=lon
   dat!0="lat"
   dat!1="lon"
   dat&lat=lat
   dat&lon=lon

  ; Subset for Southern California
   lat1 = 100
   lat2 = 500
   lon1 = 0
   lon2 = 250
   ;
   lat_out = lat(lat1:lat2)
   lon_out = lon(lon1:lon2)
   dat_out = dat(lat1:lat2,lon1:lon2)
   dat_out at long_name    = "Daily Total Precipitation"
   dat_out at units        = "mm"
   dat_out at _FillValue   = NODATA
   dat_out at missingvalue = NODATA
   nlat    = dimsizes(lat_out)
   nlon    = dimsizes(lon_out)
   ;
   ; Write out to netcdf file
   diro      = dir_out
   filo      = file_out+".nc"
   system("/bin/rm -f " + diro + filo)
   fout  = addfile (diro + filo, "c")
   setfileoption(fout,"DefineMode",True)
   fAtt               = True
   fAtt at title         = ""
   fAtt at source_file   =  ""
   fAtt at Conventions   = "None"
   fAtt at creation_date = systemfunc("date")
   fileattdef( fout, fAtt )
   dimNames = (/"time","lat","lon" /)
   dimSizes = (/ 1,nlat,nlon /)
   dimUnlim = (/ True,False,False /)
   filedimdef(fout,dimNames,dimSizes,dimUnlim)
   filevardef(fout,"lat" ,typeof(lat_out),(/"lat"/))
   filevardef(fout,"lon",typeof(lon_out),(/"lon"/))
   filevardef(fout,"ppt",typeof(dat_out),(/"lat","lon"/))
   ;
   fout->lat                 = lat_out
   fout->lon                 = lon_out
   fout->ppt                 = dat_out
   delete(fout)
   delete(lat)
   delete(lon)
   delete(dat)
   delete(dat_out)
end