ADCP Info
Plotting ADCP Data vs. Model Output
Plotting vertical profiles of currents along transects across the shelf requires both sigma coordinate interpolation and plotting the results of the regular, interplated grid. NCL offers solutions for both.
Sigma Coordinate Interpolation
CDO Example
The Climate Data Operators (CDO) at:
https://code.zmaw.de/projects/cdo
offer a nifty way to vertically interpolate our files with the cdo intlevel operator. The manual explains this at:
https://code.zmaw.de/embedded/cdo/1.5.2/cdo.html#x1-5200002.12.6.
A minimal example wherein we linearly interpolate our file H_iasroms_201001011200-time0.nc, which contains just the first time step for prototyping purposes, is:
cdo intlevel,0,20,40,60,80,100 H_iasroms_201001011200-time0.nc H_iasroms_201001011200-time0-newlev.nc
where part of the vertical field from 0 to 3500 m is linearly interpolated to 0,20,40,60 and 100 m.
NCL Example
A sigma coordinate interpolation example can be found at:
http://www.ncl.ucar.edu/Applications/sigma.shtml
with the actual code being:
;************************************
; sigma_1.ncl
;************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;************************************
begin
;************************************
; open netcdf file and read in data
; this dataset is a 2D shelf model with the
; y-axis being distance in m from the bottom
; of the shelf to the surface, and the
; x-axis being distance along the shelf
;***********************************
in = addfile("./shelf_grid.nc","r")
y = in->y ; dimensioned (dist,sigma)
d = in->DENSITY
;************************************
; create new depth array
; we will interpolate the sigma coords
; to these values
;***********************************
depth = new(120,"float")
depth(0:19) = ispan(0,99,5)
depth(20:119) = ispan(100,199,1)
;************************************
; interpolate y (sigma) to depth
; arguments (sigma values for column,data values for column,
;***********************************
nosigma = new((/4000,dimsizes(depth)/),"float")
; because each column has a unique set of coordinates, we
; must interpolate them individually. The total number of
; interpolations here (4000) slows the execution down a bit.
nosigma = linint1(y,d,False,depth,0)
nosigma!0 = "dist"
nosigma!1 = "depth"
nosigma_s = nosigma(depth|:,dist|:)
nosigma_s@long_name = d@long_name ; long_name used auto by plot template
;************************************
; create plot
;***********************************
wks = gsn_open_wks("ps","sigma") ; open a ps file
gsn_define_colormap(wks,"gui_default") ; choose colormap
res = True ; plot mods desired
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = False ; turn off contour lines
res@gsnSpreadColors = True ; use full colormap
res@lbOrientation = "vertical" ; vertical label bar
res@sfXArray = y&dist ; x axis array
res@sfYArray = depth ; y axis array
res@tiYAxisString = "Dist from bottom to surface (m)"
res@tiXAxisString = "Dist from origin (km)"
res@tiMainString = "Data now in regular coordinates"
plot = gsn_csm_contour(wks,nosigma_s,res) ; create plot
;************************************
end
Working with Model 1
The Model 1 files have been copied to:
terrebonne.tamu.edu/raid2/GOMEX/STEP1/1_IASNFS_Ko
for prototyping purposes.
Creating Subsets with NCO
Since each original file is huge (14G), subset files have been extracted to minimize memory swapping during the prototype stage. The 3-D U, V and coordinate variable fields were extracted from:
H_iasroms_201001011200.nc (14G)
via the NetCDF Operators (NCO) command:
ncks -v time,eastward_sea_water_velocity,northward_sea_water_velocity H_iasroms_201001011200.nc H_iasroms_201001011200-uv.nc
to obtain the smaller file:
H_iasroms_201001011200-uv.nc (6.3G)
Just the first time step of the entire file can be extracted via:
ncks -d time,0 H_iasroms_201001011200.nc H_iasroms_201001011200-time0.nc
The depth values with which we'll be working can be found via:
ncks -v zlev H_iasroms_201001011200.nc
and are:
zlev[0]=0 m zlev[1]=20 m zlev[2]=40 m zlev[3]=60 m zlev[4]=80 m zlev[5]=100 m zlev[6]=150 m zlev[7]=200 m zlev[8]=250 m zlev[9]=300 m zlev[10]=400 m zlev[11]=500 m zlev[12]=600 m zlev[13]=700 m zlev[14]=800 m zlev[15]=900 m zlev[16]=1000 m zlev[17]=1500 m zlev[18]=2000 m zlev[19]=2500 m zlev[20]=3000 m zlev[21]=3500 m
A list of variables within the NetCDF file can be found via:
ncdump -h H_iasroms_201001011200.nc | grep float
float longitude(longitude) ; float latitude(latitude) ; float zlev(zlev) ; float time(time) ; float sea_floor_depth_below_sea_level(latitude, longitude) ; float sea_surface_height_above_sea_level(time, latitude, longitude) ; float eastward_sea_water_velocity_at_bottom(time, latitude, longitude) ; float northward_sea_water_velocity_at_bottom(time, latitude, longitude) ; float air_pressure_at_sea_level(time, latitude, longitude) ; float eastward_sea_water_velocity(time, zlev, latitude, longitude) ; float northward_sea_water_velocity(time, zlev, latitude, longitude) ; float sea_water_potential_temperature(time, zlev, latitude, longitude) ; float sea_water_salinity(time, zlev, latitude, longitude) ;
Before starting the ncl interpreter, we must supply the NCARG_ROOT location as:
export NCARG_ROOT=/opt/NCARG
and then start the interpreter as:
ncl
which brings up the following:
Copyright (C) 1995-2011 - All Rights Reserved University Corporation for Atmospheric Research NCAR Command Language Version 6.0.0 The use of this software is governed by a License Agreement. See http://www.ncl.ucar.edu/ for more details. ncl 0>
Our script looks like this:
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
in = addfile("./H_iasroms_201001011200-uv.nc","r")
; -----------------------------------------------------------
; Create a transect and interpolate ut along that transect.
; -----------------------------------------------------------
; Use just the first time step for our prototype.
ut = in->eastward_sea_water_velocity(0,:,:,:)
; Convert any NaNs in the NetCDF file to the MissingValue so they will be
; ignored by the plotting program rather than crash it.
if (any(isnan_ieee(ut))) then
value = -9999.
replace_ieeenan(ut,value,0)
ut@_FillValue = value
end if
leftlat = 29.0
rightlat = 28.0
leftlon = -92.0
rightlon = -91.5
; western MCH grid boundary
;leftlat = 27.585
;rightlat = 29.410
;leftlon = -94.595
;rightlon = -94.644
; eastern MCH grid boundary
;leftlat = 29.103
;rightlat = 30.276
;leftlon = -87.839
;rightlon = -87.808
npts = 500
; gc_latlon - interpolates points along a great circle on the surface of the globe
;
http://www.ncl.ucar.edu/Document/Functions/Built-in/gc_latlon.shtml
; The -2 indicates distances returned in degrees, longitudes from -180 to 180.
dist = gc_latlon(leftlat,leftlon,rightlat,rightlon,npts,-2)
; linint2_points - interpolates from a grid to arbitrarily specified coordinate pairs
;
http://www.ncl.ucar.edu/Document/Functions/Built-in/linint2_points.shtml
trans = linint2_points(ut&longitude,ut&latitude,ut,False,dist@gclon,dist@gclat,2)
copy_VarAtts(ut,trans) ; copy attributes
trans!0 = "z_ut" ; create named dimension and assign
trans&z_ut = ut&zlev ; coordinate variable for 0th dimension only
;********************************
; create plot
;********************************
wks = gsn_open_wks("ps","trans")
gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; choose colormap
res = True ; plot mods desired
res@tmXBMode = "Explicit" ; explicitly label x-axis
res@tmXBValues = (/points(0),points(npts-1)/) ; points to label
; label values
res@tmXBLabels = (/leftlat +", "+leftlon,rightlat+", "+rightlon/)
res@cnFillOn = True ; turn on color
res@lbLabelAutoStride = True ; nice label bar label stride
res@gsnSpreadColors = True ; use full range of colormap
res@cnLinesOn = False ; turn off countour lines
res@lbOrientation = "vertical" ; vertical label bar
res@pmLabelBarOrthogonalPosF = -0.05 ; move label bar closer to plot
res@tiMainString = "Transect" ; add title
res@tiXAxisString = "lat/lon along transect"
res@trYReverse = True ; reverse y axis
res@trXReverse = True ; reverse x axis (neg longitudes)
res@cnLevelSpacingF = 0.05 ; set contour spacing
plot = gsn_csm_contour(wks,trans,res) ; create plot
frame(wks)
end
Transect Plotting
NCL Example A transect plotting example can be found at:
http://www.ncl.ucar.edu/Applications/transect.shtml
with the script being:
;************************************
; CSM_Graphics: transect_1.ncl
;************************************
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
in = addfile("h_avg_Y0191_D000.00.nc","r")
t = in->T(0,:,:,:)
;************************************
; calculate great circle along transect
;************************************
leftlat = -60.0
rightlat = -30.0
leftlon = -60.0
rightlon = 20.0
npts = 100 ; number of points in resulting transect
dist = gc_latlon(leftlat,leftlon,rightlat,rightlon,npts,2)
points = ispan(0,npts-1,1)*1.0
;********************************
; interpolate data to great circle
;********************************
trans = linint2_points(t&lon_t,t&lat_t,t,True,dist@gclon,dist@gclat,2)
copy_VarAtts(t,trans) ; copy attributes
trans!0 = "z_t" ; create named dimension and assign
trans&z_t = t&z_t ; coordinate variable for 0th dimension only
;********************************
; create plot
;********************************
wks = gsn_open_wks("ps","trans") ; open a ps file
gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; choose colormap
res = True ; plot mods desired
res@tmXBMode = "Explicit" ; explicitly label x-axis
res@tmXBValues = (/points(0),points(npts-1)/) ; points to label
; label values
res@tmXBLabels = (/leftlat +", "+leftlon,rightlat+", "+rightlon/)
res@cnFillOn = True ; turn on color
res@lbLabelAutoStride = True ; nice label bar label stride
res@gsnSpreadColors = True ; use full range of colormap
res@cnLinesOn = False ; turn off countour lines
res@lbOrientation = "vertical" ; vertical label bar
res@pmLabelBarOrthogonalPosF = -0.05 ; move label bar closer to plot
res@tiMainString = "Transect" ; add title
res@tiXAxisString = "lat/lon along transect"
res@trYReverse = True ; reverse y axis
; res@trXReverse = True ; reverse x axis (neg longitudes)
res@cnLevelSpacingF = 1.0 ; set contour spacing
plot = gsn_csm_contour(wks,trans,res) ; create plot
;********************************
; show transect on a map
;********************************
i = NhlNewColor(wks,0.8,0.8,0.8) ;add gray to continents
mres = True ; plot mods desired
mres@gsnFrame = False ; don't turn page yet
mres@gsnDraw = False ; don't draw yet
mres@tiMainString= "Transect Location" ; title
map = gsn_csm_map_ce(wks,mres) ; create map
; add polyline
pres = True ; polyline mods desired
pres@gsLineColor = "red" ; color of lines
pres@gsLineThicknessF = 2.0 ; line thickness
gsn_polyline(wks,map,(/leftlon,rightlon/),(/leftlat,rightlat/),pres)
draw(map)
frame(wks)
end