How to warp an image across a portion of a map projection

I have an image I would like to warp to a map projection in IDL. How do I do this? 

Specifically, I would like to have control over the output size of the resulting window. I have tried MAP_SET and MAP_IMAGE, and the results seem difficult to control.

First, consider the following:

MAP_IMAGE is a function which returns an image that is warped to fit a current map projection. Therefore, MAP_SET must be called before MAP_IMAGE in order to set up that projection.

Also consider that MAP_SET takes care of the data to device coordinate conversion. MAP_SET also establishes the projection types expressed in the bounds of the customized projection you choose.

KEY TIP:
A common problem seen is that the boundries of the image to be mapped, and the resulting projection veiw seem difficult to control. The most direct way to have control of these boundries is to set the LIMIT keyword to BOTH MAP_SET and to MAP_IMAGE. See this in the example below.

In MAP_SET, this LIMIT is a vector that specifies the boundries to be mapped.

In MAP_IMAGE, this LIMIT represents the boundries im map coordinates of the warped image (based in terms of the first and last rows and columns in the image)

If the LIMIT is not set in MAP_IMAGE, the default behavior that is seen is that the warped image might be displayed outside the bounds that are set to the limit in MAP_SET.

Below is the sequence of events that should occur in order to warp your image to the map projection you are establishing. First, in English, and then in IDL code.

1) The procedure is called image_warp2_map

2) Set visual class to use decomposed color if not in pseudo color. Also set RETAIN = 2 to make IDL handle backing store.

3) Make some image data.

4) Set the boundries for the LIMIT vector
Minlat=-65
Maxlat=65
Minlon = 160     ;Left edge is 160 East
Maxlon = -70 + 360    ;Right edge is 70 West = +290.

5) Define the limit1 vector to be:
limit1 = [minlat, minlon, maxlat, maxlon]

6) Make some more data, and define loncenter
loncenter= -140
idata = replicate(100, 50, 50)

7) Use the WINDOW procedure to explicitly define the XSIZE and YSIZE of the display window. window, 1, xsize = 255, ysize = 255

8) Call MAP_SET. In the example, you will see that: The LAT to be mapped to the center of the image is at 0 degrees. The LON to be mapped to the center is the predefined loncenter from step 6 above. The default angle of rotation is 0 degrees. The map projection is Orthographic, and the boundaries of the device window are set to limit1.

Note here that MAP_SET will figure out the conversion to make latmin to latmax go from 0-255 in device coordinates, and apply this conversion to the image data.

MAP_SET, 0, loncenter, 0, /ortho, limit=limit1

9) Call MAP_IMAGE. The resulting warped image will start at startx and starty. Note that these starting points are off from the lower left corner of the display window by the amount (xstart, ystart). These values are returned to MAP_IMAGE in device coordinates, therefore keeping the image from touching the lower horizontal and left vertical boundries of the image window.

XSIZE and YSIZE were defined in the call to WINDOW. Note here that the image will scale to this window size. Therefore, changing the size in the call to the WINDOW procedure will give you control of the output size of the warped image.

Then, be sure to set the lat and lon max and min to appropriate subscript of the limit1 vector. This will ensure that the boundries of the warped image fit in the inside of the display window and also within the boundries of the map projection.

result=MAP_IMAGE(idata,startx,starty,xsize,ysize,$
latmin=limit1(0),latmax=limit1(2),$
lonmin=limit1(1),lonmax=limit1(3))

10) Next, you can display the warped image on the map projection at the proper position (which starts at (startx, starty) rather than at (0,0))

11) Lastly, draw the continents and the grid lines over the map.

Example:

 

PRO make_a_map2 
DEVICE, DECOMPOSE = 0, RETAIN = 2
LOADCT,5 
IMAGE = BYTSCL(SIN(DIST(40)/10)) Minlat = -65
Maxlat = 65
Minlon = 160       ;Left edge is 160 East
Maxlon = -70 + 360   ;Right edge is 70 West = +290. 
limit1 = [minlat, minlon, maxlat, maxlon] 
loncenter= -140
idata = REPLICATE(100, 50, 50) 
WINDOW, 1, XSIZE = 255, YSIZE = 255  
;uncomment the line below to see a large window;
with the map scaled to the window accordingly 
;WINDOW, 1, XSIZE = 500, YSIZE = 500 
;uncomment the line below to see a small window
;with the map scaled to the window accordingly 
;WINDOW, 1, XSIZE = 200, YSIZE = 200 
MAP_SET, 0, loncenter, 0, /ORTHO, LIMIT=limit1 
result=MAP_IMAGE(idata,startx,starty,xsize,ysize,$     
 latmin=limit1(0),latmax=limit1(2),$    
 lonmin=limit1(1),lonmax=limit1(3)) 
tv,result,startx,starty   
MAP_CONTINENTS, /COASTSMAP_GRID,latdel=10,londel=10, /LABEL,     $    
 lonlab=limit1(0)+.5, latlab=limit1(1),$    
 glinestyle=0    
END

Review on 12/31/2013 MM