[ncl-talk] Reading PRISM 4km data onto regular lat/lon grid
Dave Allured - NOAA Affiliate
dave.allured at noaa.gov
Thu Oct 1 11:51:51 MDT 2015
Scott,
Don't use fspan for this purpose. The following method guarantees precise
increments to the nearest possible values in double precision. This method
requires that the intended theoretical increments are exact fractions of
one degree, which I believe is the correct assumption for 4 km PRISM.
opt_integer = 3 ; get exact integer divisors
lats_per_degree = round (1.0 / YDIM, opt_integer)
lons_per_degree = round (1.0 / XDIM, opt_integer)
dlats = ULYMAP - todouble (ispan (0, nlats-1, 1)) / lats_per_degree
dlons = ULXMAP + todouble (ispan (0, nlons-1, 1)) / lons_per_degree
Once you have the nearest possible values in double precision, you can then
convert to single precision values that are accurate within one bit, if you
like.
lats = tofloat (dlats)
lons = tofloat (dlons)
--Dave
On Thu, Oct 1, 2015 at 10:53 AM, Scott Capps <scott at allvertum.com> wrote:
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151001/86410986/attachment.html
More information about the ncl-talk
mailing list