PRO SEAMASK, Position=position

   ; POSITION: Four-element array indicating position of map in the output window.
   IF N_Elements(position) EQ 0 THEN position=[0.0,0.0,1.0,1.0]; Fill up window.

      ; Program needs non-decomposed color.

   IF (!D.Flags AND 256) NE 0 THEN Device, Decomposed=0, Get_Decomposed=theState

      ; Pick a seed, so you see what I see. Create data.

   seed = -5L
   lat = RANDOMU(seed, 50) * (24) + 23
   lon = RANDOMU(seed, 50) * 50.0 - 123
   data = RANDOMU(seed, 50) * 100

      ; Colors for plot.

   landColor    = !D.Table_Size-2
   contourColor = !D.Table_Size-3
   seaColor     = !D.Table_Size-4
   dataColor    = !D.Table_Size-5
   usaColor     = !D.Table_Size-6
   contColor    = !D.Table_Size-7
   LOADCT, 0
   TVLCT, FSC_Color('Wheat', /Triple), landColor
   TVLCT, FSC_Color('Firebrick', /Triple), contourColor
   TVLCT, FSC_Color('Sky Blue', /Triple), seaColor
   TVLCT, FSC_Color('Coral', /Triple), dataColor
   TVLCT, FSC_Color('Dark Gray', /Triple), usaColor
   TVLCT, FSC_Color('Indian Red', /Triple), contColor

      ; The TVRD's require an 8-bit image. This is most
      ; easily done in the Z-graphics buffer.

   thisDevice = !D.Name
   Set_Plot, 'Z', /COPY
   Device, Set_Resolution=[500, 350]
   Erase, Color=0

      ; Set up the map projection of the US continent.

   mapPosition = position
   MAP_SET, 15, -87, 0, LIMIT=[23,-125,50,-67], /NOBORDER, $
      /MERCATOR, /NOERASE, POSITION=mapPosition, $
      XMARGIN=0, YMARGIN=0

      ; Draw the filled continental outlines.

   MAP_CONTINENTS, /Fill_Continents, Color=255

      ; Use a snapshot of this image for a sea mask.

   snap = TVRD( mapPosition[0]*!D.X_Size, mapPosition[1]*!D.Y_Size, $
      (mapPosition[2]-mapPosition[0])*!D.X_Size, $
      (mapPosition[3]-mapPosition[1])*!D.Y_Size)
   mask = (snap EQ 0)
   sea =  mask * seaColor

      ; Draw the continental outlines and the USA state boundaries.

   TV, sea, mapPosition[0], mapPosition[1], /Normal
   MAP_CONTINENTS, /Fill, Color=landColor
   MAP_CONTINENTS,  Color=contColor
   MAP_CONTINENTS, /USA, Color=usaColor

      ; Plot the random data locations.

   PLOTS, lon, lat, PSYM=4, COLOR=dataColor

      ; Grid the irregularly spaced data.

   gridData= SPH_SCAT(lon, lat, data, $
      BOUNDS=[-125., 23., -67., 50.], GS=[0.5,0.5], BOUT=bout)

      ; Calculate xlon and ylat vectors corresponding to gridded data.

   s = SIZE(gridData)
   xlon = FINDGEN(s(1))*((bout(2) - bout(0))/(s(1)-1)) + bout(0)
   ylat = FINDGEN(s(2))*((bout(3) - bout(1))/(s(2)-1)) + bout(1)

      ; Put the contours on the map.

   CONTOUR, gridData, xlon, ylat, /OVERPLOT, NLEVELS=12, C_COLOR=contourColor, $
      C_LABEL=Replicate(1, 12), /DOWNHILL

      ; Take a snapshot of what is currently in the window.

   contourImage = TVRD( mapPosition[0]*!D.X_Size, mapPosition[1]*!D.Y_Size, $
      (mapPosition[2]-mapPosition[0])*!D.X_Size, $
      (mapPosition[3]-mapPosition[1])*!D.Y_Size)

      ; Replace the sea pixels in the contour image with pixels from the sea image.

   theSea = Where(sea GT 0, count)
   finalImage = contourImage
   IF count GT 0 THEN finalImage[theSea] = sea[theSea]

      ; Back to the display device and display the two images.

   Set_Plot, thisDevice
   IF (!D.Flags AND 256) NE 0 THEN Window, XSIZE=500, YSIZE=350, Title='Non-Masked Map'
   TV, contourImage, mapPosition[0], mapPosition[1], /Normal
   IF (!D.Flags AND 256) NE 0 THEN Window, 1, XSIZE=500, YSIZE=350, Title='Sea-Masked Map'
   TV, finalImage, mapPosition[0], mapPosition[1], /Normal

      ; Set the decomposed state back.

   IF N_Elements(theState) NE 0 THEN Device, Decomposed=theState
   END