Sea Masks on Map Projections
QUESTION: I have some data that I want to display as contours on a map projection of a continent. But the contours extend over the water and the data pertains only to the land mass. What I would like to do is mask all the contours overlaying water out of my display. But I can't find an easy way to do this. Can it be done in IDL?
ANSWER: If by "easy way" you mean by setting a keyword or calling an IDL command to do it, then no, there is no easy way. But it can be done, however. :-)
The essential technique is to create a sea mask. You might think this would be done with a routine like MapContinents and a WaterOnly keyword, but life is not so simple. We do need to use MapContinents, however, because it is possible to create a land mask with that command. The sea mask is, of course, the inverse of the land mask.
There are probably a number of different ways to do this. I'm going to show you a method that involves the Z-graphics buffer, simply because this is easy and straightforward. I want the masks to be simple 8-bit images, and these are easy to create in an 8-bit device like the Z-graphics buffer, and more difficult to create on a 24-bit device. (Although it certainly can be done.)
I have demonstrated the technique in an example program, named Seamask, which I will take you through in detail here. The program takes some simulated random data at random locations over the United States land mass, grids the data so it can be contoured, and displays it on a map projection. You will need to download the FSC_Color program from the Coyote Library to be able to run the example program
The first part of the program simply creates the random data that will be contoured and sets up the colors to be used. It looks like this.
PRO SEAMASK
; 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) * 65.0 - 120
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
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
I want the map to be an 8-bit image. This is most easily done in the Z-graphics buffer, which is an 8-bit device, where the TVRD command will return an 8-bit image. This can also be done on a 24-bit device, but you will have to be a bit more careful. I enter the Z-graphics buffer and get it set up with these commands.
; 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
Next, I set up my map projection space. This is most easily done by ensuring the map projection extends to the edges of the display window, but it doesn't have to do this. If you prefer not to do this, you will probably have to use the POSITION keyword to position your map projection and then calculate exactly where the image is in the larger window so that it can be extracted for the mask. I've made it easy on myself by using the entire window for the mask, which I can capture using TVRD without any arguments.
; Set up the map projection of the US continent.
MAP_SET, 15, -87, 0, LIMIT=[23,-125,50,-67], /NOBORDER, $
/MERCATOR, /NOERASE, POSITION=[0,0,1,1], $
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()
mask = (snap EQ 0)
sea = mask * seaColor
Now I draw the essential elements of my map projection, including the sea, land, state boundaries, and the data and contour locations themselves.
; Draw the sea, the continental outlines, and the USA state boundaries.
TV, sea
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
I want to save what is currently in the window as the "without sea mask" result.
; Take a snapshot of what is currently in the window.
contourImage = TVRD()
What I need to do now is replace the sea pixels in the contour image with the corresponding pixels from the sea image. The code looks like this.
; 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]
Finally, we return from the Z-graphics buffer and display the two images.
; 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
IF (!D.Flags AND 256) NE 0 THEN Window, 1, XSIZE=500, YSIZE=350, Title='Sea-Masked Map'
TV, finalImage
; Set the decomposed state back.
IF N_Elements(theState) NE 0 THEN Device, Decomposed=theState
END
You can see the results of the Seamask program in the two figures below.
Copyright © 2003 David W. Fanning
Last Updated 16 July 2003