Masking Question

Hello Again:

Thanks for the help you have given me in the past.

This question is about masking. I've included a script (see below). It
works, but I don't know how to proceed from here. The goal is:

1) Give all land one color.
2) Give each ocean point a color depending on the sea surface
temperature of that point.
3) Mask out all cloudy points over the ocean.

The script does all this. So far so good. From here I'm lost. I now
want to:

4) Give all the ocean--except the points colored by a sea surface
temperature--one color (example: gray).

Will I have to restructure the masks? I'm not sure. Your help is much
appreciated.

---A Million Thanks
---Tony Conrad
---NSOF/NOAA
---Suitland, Maryland, USA

;*********************************
; alaskapp_4.ncl
; This uses the 2D lat/lon array plotting procedure available in
; NCL version 4.2.0.a028 and higher (version 4.2.0.a025 contains
; the 2D plotting technique, but a bug that prevented this from
; being used with EOSDIS data was not corrected until version 4.2.0.a028).
;*********************************
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
;*********************************
; read in HDF file
;*********************************
  f1 =
addfile("/usr/data/modisRedo/223-Day/A2008223162500.L2_LAC_SST.hdf", "r")
; f1 =
addfile("/usr/data/modisRedo/223-Day/A2008223000000.L2_LAC_SST.hdf", "r")
;*************************************************************
; read in sst data
;*************************************************************
  ch1_int = f1->sst
; printVarSummary(ch1_int)
  ch1 = ch1_int*ch1_int_at_slope(0) ; orig is integer, make float
; printVarSummary(ch1)
; printMinMax(ch1, True)
;*********************************
; read in second HDF file
;*********************************
  f2 =
addfile("/usr/data/modisRedo/223-Day/MYD03.A2008223.1625.005.2008224143628.hdf",
"r")
; f2 =
addfile("/usr/data/modisRedo/223-Day/MYD03.A2008223.0000.005.2008223165522.hdf",
"r")
;*************************************************************
; read in navigation data about the sst data
;*************************************************************
  lat2d = f2->Latitude
  lon2d = f2->Longitude

; printVarSummary(lat2d)
; printMinMax(lat2d, True)
; printVarSummary(lon2d)
; printMinMax(lon2d, True)
;*********************************
; read in third HDF file
;*********************************
  f3 =
addfile("/usr/data/modisRedo/223-Day/MYD35_L2.A2008223.1625.005.2008224153726.hdf",
"r")
;*************************************************************
; read in cloud mask data about the sst data
;*************************************************************
  cloud = f3->Cloud_Mask(0,:,:)
; printVarSummary(cloud)
; print(cloud)
;*************************************************************
; isolate least significant bits (cloud mask) in cloud mask data
;*************************************************************
  c_bits = dim_gbits( cloud, 5, 3, 5, 1354 ) ; isolate three bits
; printVarSummary(c_bits)
; print(c_bits)
;*************************************************************
; isolate most significant bits (water/land bits) in cloud mask data
;*************************************************************
  w_bits = dim_gbits( cloud, 0, 2, 6, 1354 ) ; isolate two bits
; printVarSummary(w_bits)
; print(w_bits)
;*************************************************************
; mask out cloudy sst data and sst data over land
;*************************************************************
; ch2 = mask( ch1, (c_bits.ne.1 .and. ch1.gt.-150.00), True )
  ch2 = mask( ch1, (c_bits.ne.1 .and. w_bits.eq.0), True )

  ch2_at_lat2d = lat2d ; assign nav data to sst data
  ch2_at_lon2d = lon2d
;*************************************************************
; assign named dimensions
;*************************************************************
  ch2_at_long_name = "Sea Surface Temp"
  ch2_at_units = "degrees Celsius"
;*************************************************************
; create plot
;*************************************************************
  wks = gsn_open_wks("x11", "alaska")
; gsn_define_colormap(wks, "ncview_default") ; choose colormap
  gsn_define_colormap(wks, "rainbow+white+gray") ; choose colormap
 
  res = True
  res_at_cnFillOn = True ; color Fill
  res_at_cnFillMode = "RasterFill" ; Raster Mode
  res_at_cnLinesOn = False ; Turn off contour lines
 
  res_at_cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  res_at_cnMinLevelValF = -50.0 ; set min contour level
  res_at_cnMaxLevelValF = 30.0 ; set max contour level
  res_at_cnLevelSpacingF = 1.0 ; set contour spacing
 
  res_at_gsnSpreadColors = True ; use full colormap
  res_at_gsnSpreadColorStart = 14 ; start at color 14
  res_at_gsnSpreadColorEnd = 237 ; end at color 237
  res_at_lbLabelStride = 10 ; every 10nth label bar label
  res_at_lbOrientation = "Vertical" ; vertical label bar
 
  res_at_mpProjection = "Stereographic" ; polar stereographic
  res_at_mpLimitMode = "LatLon" ; map bounded by lat/long
  res_at_mpMinLatF = min(lat2d) ; choose lat/lon map limits
  res_at_mpMaxLatF = max(lat2d)
  res_at_mpMinLonF = min(lon2d)
  res_at_mpMaxLonF = max(lon2d)
  res_at_gsnMaximize = True ; maximize (without
                                            ; overflowing) the plot area

;*************************************************************
; turn on lat/lon labeling
;*************************************************************
  res_at_pmTickMarkDisplayMode = "Always" ; turn on tickmarks
  res_at_mpGridAndLimbOn = True ; turn on lat/lon lines
  res_at_mpGridLatSpacingF = 5 ; spacing for lat lines
  res_at_mpGridLonSpacingF = 10 ; spacing for lon lines

  res_at_mpFillDrawOrder = "PostDraw" ; draw land last
  res_at_mpLandFillColor = 191 ; color for land
  res_at_mpDataBaseVersion = "Ncarg4_1" ; higher res coastline
  res_at_gsnAddCyclic = False ; needed b/c data not global
 
  res_at_tiMainString = "Modis data on a polar stereographic projection"

 
  plot = gsn_csm_contour_map(wks, ch2, res) ; create plot
 
end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Sep 09 2008 - 09:06:15 MDT