Re: problems with writing out netcdf files and plotting with lat/lon2d

Kerrie,

NetCDF files do not recognize 2D attributes. When you write "topo" to a NetCDF file with the attached lat2d/lon2d attributes, they get converted to 1D arrays.

Hence, when you try to read it back in, topo@lat2d and topo@lon2d are now 1D arrays, and NCL gives you the error that your lat/lon coordinate arrays don't match the size of your data array.

Using "lat2d" and "lon2d" as special attributes to "topo" are for plotting only. They should not really be used for writing data back out to a file.

I actually prefer *not* using the "topo@lat2d"/"topo@lon2d" method, and instead doing something like this, where you use the special sfXArray and sfYArray resources for setting lat2d and lon2d.

Note that in the code below, I write out lon2d and lat2d separately to the file, since they are no longer attributes of "topo".

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/wrf/WRFUserARW.ncl"

begin
                            
infile=addfile("/Network/p1/WRF_IPCC_OUTPUT/OUTPUT/wrfout_d01_1998-01-03_00:00:00.nc","r")
XLAT=wrf_user_getvar(infile,"XLAT",-1)
XLONG=wrf_user_getvar(infile,"XLONG",-1)
HGT=wrf_user_getvar(infile,"HGT",-1)

zoomlats=(/15.,35./)
zoomlons=(/-120.,-80./)
loc1=wrf_user_ll_to_ij(infile,zoomlons,zoomlats,True)
  x_start = loc1(0,0) - 1
  x_end = loc1(0,1) - 1
  y_start = loc1(1,0) - 1
  y_end = loc1(1,1) - 1

lat2d=XLAT(0,y_start:y_end,x_start:x_end)
lon2d=XLONG(0,y_start:y_end,x_start:x_end)

topo=HGT(0,y_start:y_end,x_start:x_end)
delete_VarAtts(topo,-1)
topo@units="m"
;;;topo@lat2d=lat2d
;;;topo@lon2d=lon2d

wks = gsn_open_wks("x11","test")

res=True
res@mpMinLatF=min(topo@lat2d)
res@mpMaxLatF=max(topo@lat2d)
res@mpMinLonF=min(topo@lon2d)
res@mpMaxLonF=max(topo@lon2d)

;---Define lat/lon coordinates of topo
res@sfXArray = lon2d
res@sfYArray = lat2d
res@gsnAddCyclic = False ; necessary for regional data

plot1=gsn_csm_contour_map(wks,topo,res) ;### THIS PLOT IS CORRECT

system("rm -f /Volumes/ocean6/WRF_HadCM3_downscale/WRF_HadCM3.TOPO.nc")
cdf_file=addfile("/Volumes/ocean6/WRF_HadCM3_downscale/WRF_HadCM3.TOPO.nc","c")
    cdf_file->topo=topo
    cdf_file->lat2d = lat2d
    cdf_file->lon2d = lon2d

infile=addfile("/Volumes/ocean6/WRF_HadCM3_downscale/WRF_HadCM3.TOPO.nc","r")
; topo=infile->topo
    topo=wrf_user_getvar(infile,"topo",-1)

plot2=gsn_csm_contour_map(wks,topo,res) ;### THIS PLOT IS WRONG

end

On Oct 17, 2011, at 10:58 PM, Kerrie Geil wrote:

> Hi,
>
> I can read in topography directly from a WRF output file and it will plot correctly with lat2d and lon2d
> I then subset the topo to a smaller coverage area and write the subsetted topo to a new netcdf file.
> When I read the topo from the new netcdf file and try to plot, the plot is incorrect and I get these errors for both lat and lon:
>
> (0) is_valid_latlon2d_attr: Warning: The 'lat2d' attribute must either be the same dimension sizes as the data, or one element larger in both directions. Your data will most likely not be overlaid on the map correctly.
> (0) check_for_y_lat_coord: Warning: Data either does not contain a valid latitude coordinate array or doesn't contain one at all.
> (0) A valid latitude coordinate array should have a 'units' attribute equal to one of the following values:
> (0) 'degrees_north' 'degrees-north' 'degree_north' 'degrees north' 'degrees_N' 'Degrees_north' 'degree_N' 'degreeN' 'degreesN' 'deg north'
>
> Is there something else I need to write to the new netcdf file for the second plot to work correctly? Thanks for the help.
>
> My code is:
>
> 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/wrf/WRFUserARW.ncl"
>
>
> begin
>
> infile=addfile("/Network/p1/WRF_IPCC_OUTPUT/OUTPUT/wrfout_d01_1998-01-03_00:00:00.nc","r")
> XLAT=wrf_user_getvar(infile,"XLAT",-1)
> XLONG=wrf_user_getvar(infile,"XLONG",-1)
> HGT=wrf_user_getvar(infile,"HGT",-1)
>
> zoomlats=(/15.,35./)
> zoomlons=(/-120.,-80./)
> loc1=wrf_user_ll_to_ij(infile,zoomlons,zoomlats,True)
> x_start = loc1(0,0) - 1
> x_end = loc1(0,1) - 1
> y_start = loc1(1,0) - 1
> y_end = loc1(1,1) - 1
>
> lat2d=XLAT(0,y_start:y_end,x_start:x_end)
> lon2d=XLONG(0,y_start:y_end,x_start:x_end)
>
> topo=HGT(0,y_start:y_end,x_start:x_end)
> delete_VarAtts(topo,-1)
> topo@units="m"
> topo@lat2d=lat2d
> topo@lon2d=lon2d
>
> wks = gsn_open_wks("x11","test")
>
> res=True
> res@mpMinLatF=min(topo@lat2d)
> res@mpMaxLatF=max(topo@lat2d)
> res@mpMinLonF=min(topo@lon2d)
> res@mpMaxLonF=max(topo@lon2d)
>
> plot1=gsn_csm_contour_map(wks,topo,res) ;### THIS PLOT IS CORRECT
>
>
> system("rm -f /Volumes/ocean6/WRF_HadCM3_downscale/WRF_HadCM3.TOPO.nc")
> cdf_file=addfile("/Volumes/ocean6/WRF_HadCM3_downscale/WRF_HadCM3.TOPO.nc","c")
> cdf_file->topo=topo
>
> infile=addfile("/Volumes/ocean6/WRF_HadCM3_downscale/WRF_HadCM3.TOPO.nc","r")
> ; topo=infile->topo
> topo=wrf_user_getvar(infile,"topo",-1)
>
> plot2=gsn_csm_contour_map(wks,topo,res) ;### THIS PLOT IS WRONG
>
>
> end
>
> --
> Kerrie Geil
>
> Master's Student
> Department of Atmospheric Sciences
> University of Arizona
> PAS Building Rm 526
> 1118 E 4th Street
> Tucson, AZ 85721
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Oct 18 09:27:49 2011