<div dir="ltr">Scott,<div><br></div><div>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.</div><div><div><br></div><div><div> opt_integer = 3 ; get exact integer divisors</div><div> lats_per_degree = round (1.0 / YDIM, opt_integer)</div><div> lons_per_degree = round (1.0 / XDIM, opt_integer)</div><div><br></div><div> dlats = ULYMAP - todouble (ispan (0, nlats-1, 1)) / lats_per_degree</div><div> dlons = ULXMAP + todouble (ispan (0, nlons-1, 1)) / lons_per_degree</div></div><div><br></div><div>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.</div><div><div><br></div><div> lats = tofloat (dlats)</div><div> lons = tofloat (dlons)</div></div><div><br></div><div>--Dave</div><div class="gmail_extra"><br></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Oct 1, 2015 at 10:53 AM, Scott Capps <span dir="ltr"><<a href="mailto:scott@allvertum.com" target="_blank">scott@allvertum.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Greetings,<br>
<br>
I am trying to read in PRISM 4km data in bil format and output on the<br>
regular lat/lon grid, saving as a NetCDF file. I eventually want to<br>
mask the data using a shapefile and want to keep it on a regular lat/lon<br>
grid.<br>
<br>
The script works fine but I am concerned about the issues below where<br>
the dx and dy do not equate to 0.04166667 degrees (See lines enclosed<br>
by: *************). I can provide the PRISM data when needed. Any help<br>
is appreciated.<br>
<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"<br>
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"<br>
; USAGE<br>
; ncl 'dir_in="./"' 'fil_in="PRISM_ppt_stable_4kmD2_19820210_bil.bil"'<br>
'dir_out="./"' 'file_out="prism_socalif_ppt_4km_02101982"'<br>
prism_to_netcdf.ncl<br>
;<br>
begin<br>
;----------------------------------------------------------------------------<br>
; FROM THE HEADER FILE:<br>
;<br>
; BYTEORDER I<br>
; LAYOUT BIL<br>
; NROWS 621<br>
; NCOLS 1405<br>
; NBANDS 1<br>
; NBITS 32<br>
; BANDROWBYTES 5620<br>
; TOTALROWBYTES 56201<br>
; PIXELTYPE FLOAT<br>
; ULXMAP -125<br>
; ULYMAP 49.9166666666687<br>
; XDIM 0.04166666666667<br>
; YDIM 0.04166666666667<br>
; NODATA -9999<br>
;----------------------------------------------------------------------------<br>
;----------------------------------------------------------------------------<br>
;----------------------------------------------------------------------------<br>
NBANDS = 1<br>
NROWS = 621<br>
NCOLS = 1405<br>
ORDER = "I"<br>
TYPE = "float"<br>
NODATA = -9999.0<br>
<br>
ULXMAP = -125.0d0<br>
ULYMAP = 49.9166666666687d0<br>
<br>
XDIM = 0.04166666666667d0<br>
YDIM = 0.04166666666667d0<br>
<br>
diri = dir_in<br>
fili = fil_in<br>
<br>
if (ORDER.ne."I") then<br>
setfileoption("bin","ReadByteOrder","LittleEndian") ; perform<br>
byte swapping<br>
end if<br>
<br>
pltType = "ps"<br>
pltDir = "./"<br>
;<br>
; I am on a littleendian machine<br>
;if (isbigendian() ) then<br>
; print("bigendian")<br>
;else<br>
; print("littleendian")<br>
;end if<br>
;<br>
do NB=0,NBANDS-1<br>
dat = fbindirread(diri+fili,NB, (/NROWS,NCOLS/), TYPE)<br>
dat@_FillValue = NODATA<br>
dat@long_name = "Daily Total Precipitation"<br>
dat@units = "mm"<br>
dat = dat(::-1,:)<br>
end do<br>
; 4km resolution (0.04166667x0.04166667 decimal degrees)<br>
lat=fspan(24.0625d0,49.9375d0,NROWS) ; Bounded by 24.0625 to 49.9375<br>
Degrees North<br>
lon=fspan(-125.0208333d0,-66.4791667d0,NCOLS) ; Bounded by<br>
-125.0208333 to -66.4791667 Degrees East<br>
; ***********************************<br>
; Should be intervals of 0.04166667 (or, close to that at least)<br>
; instead it has unequal intervals on the order of: 0.04173387096774661<br>
; ***********************************<br>
print(lat(1:NROWS-1)-lat(0:NROWS-2))<br>
; ***********************************<br>
; Should be intervals of 0.04166667<br>
; instead it has unequal intervals on the order of: 0.0416963437<br>
; ***********************************<br>
print(lon(1:NCOLS-1)-lon(0:NCOLS-2))<br>
;<br>
lat!0="lat"<br>
lat@units="degrees_north"<br>
lat&lat=lat<br>
lon!0="lon"<br>
lon@units="degrees_east"<br>
lon&lon=lon<br>
dat!0="lat"<br>
dat!1="lon"<br>
dat&lat=lat<br>
dat&lon=lon<br>
<br>
; Subset for Southern California<br>
lat1 = 100<br>
lat2 = 500<br>
lon1 = 0<br>
lon2 = 250<br>
;<br>
lat_out = lat(lat1:lat2)<br>
lon_out = lon(lon1:lon2)<br>
dat_out = dat(lat1:lat2,lon1:lon2)<br>
dat_out@long_name = "Daily Total Precipitation"<br>
dat_out@units = "mm"<br>
dat_out@_FillValue = NODATA<br>
dat_out@missingvalue = NODATA<br>
nlat = dimsizes(lat_out)<br>
nlon = dimsizes(lon_out)<br>
;<br>
; Write out to netcdf file<br>
diro = dir_out<br>
filo = file_out+".nc"<br>
system("/bin/rm -f " + diro + filo)<br>
fout = addfile (diro + filo, "c")<br>
setfileoption(fout,"DefineMode",True)<br>
fAtt = True<br>
fAtt@title = ""<br>
fAtt@source_file = ""<br>
fAtt@Conventions = "None"<br>
fAtt@creation_date = systemfunc("date")<br>
fileattdef( fout, fAtt )<br>
dimNames = (/"time","lat","lon" /)<br>
dimSizes = (/ 1,nlat,nlon /)<br>
dimUnlim = (/ True,False,False /)<br>
filedimdef(fout,dimNames,dimSizes,dimUnlim)<br>
filevardef(fout,"lat" ,typeof(lat_out),(/"lat"/))<br>
filevardef(fout,"lon",typeof(lon_out),(/"lon"/))<br>
filevardef(fout,"ppt",typeof(dat_out),(/"lat","lon"/))<br>
;<br>
fout->lat = lat_out<br>
fout->lon = lon_out<br>
fout->ppt = dat_out<br>
delete(fout)<br>
delete(lat)<br>
delete(lon)<br>
delete(dat)<br>
delete(dat_out)<br>
end<br></blockquote></div></div></div></div>