function VertexMap, Vertex, VMap, Colors,	$
	polygon=Polygon, scale=Scale, nearestneighbour=NearestNeighbour
;+
; NAME:
;	VertexMap
; PURPOSE:
;	Map a given synoptic map to a color index range.
; CATEGORY:
; CALLING SEQUENCE:
;	ColorIndices = VertexMap(Vertex, VMap, Colors,	$
;		polygon=Polygon, scale=Scale, nearestneighbour=NearestNeighbour)
; INPUTS:
;	Vertex		float array[3,n]
;						x,y,z coordinates of the vertices of a spherical
;						polygon surface (as returned by VertexSphere)
;	VMap		array[nLng,nLat]
;						synoptic map to be mapped to sphere, covering range
;						[0,360] degrees in longitude, and [-90,90] deg in latitude
; OPTIONAL INPUT PARAMETERS:
;	Colors		int array[2]
;						defines a range of color indices. VMap is rescaled
;						into the specified color range: the minimum is mapped
;						to Colors[0], the maximum to Colors[1].
;						If VMap already is an array of color indices than
;						this argument can be omitted.
;	Scale		int scalar
;						use only Scale evenly space colors between
;						Color[0] and Color[1]
;	/nearestneighbour	if set, the map value at each vertex is calculated by
;						taking the value of the neares point of VMap.
;	polygon=Polygon
;				long array[m]
;						polygon list for the spherical surface (as returned by
;						VertexSphere). See IDL Help for shade_volume for
;						description of the structure of Polygon
;
; OUTPUTS:
;	ColorIndices
;				No polygon array specified:
;				long array[n]		color indices; one for each vertex.
;									To be used in the 'shade' keyword to polyshade
;				Polygon array specied:
;				long array[m]		color indices; one for each polygon
;									To be used in the 'shade_polygon' keyword to polyshade

; OPTIONAL OUTPUT PARAMETERS:
; CALLS:
;	Interpolate, Scale2Levels
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; >	The synoptic map VMap is 2D array giving function values on a spherical surface
;	on the regular grid in longitude ([0,360] degrees) and latitude ([-90,90] degrees).
; >	The Vertex array and Polygon array describe the surface of the sphere in terms
;	of polygons.
; >	Function values at the vertices are calculated by interpolation or nearest neighbour
;	method, and converted to integer color indices, using the round function.
; MODIFICATION HISTORY:
;-
on_error, 2

Map = float(VMap)
MapMin = min(Map,max=MapMax)
if n_elements(Colors) ge 2 then		$	; Scale map to color index range specified in Colors
	if Colors[0] ne Colors[1] then	$
		Map = Colors[0]+(Colors[1]-Colors[0])*(Map-MapMin)/(MapMax-MapMin)
MapMin = min(Map,max=MapMax)
									; Convert vertex Cartesian coordinates to spherical coordinates
L = cv_coord(from_rect=Vertex, /to_sphere, /degrees)
L = L[0:1,*]							; Extract longitude and latitude

S = size(Map)
L[0,*] = L[0,*]+360.*(L[0,*] lt 0.)		; Scale longitude to [0,360]
L[0,*] = L[0,*]/360.*(S[1]-1)			; Scale longitude to [0,nLng-1]
L[1,*] = (1+L[1,*]/90.)*0.5*(S[2]-1)	; Scale latitude  to [0,nLat-1]

; Calculate values at vertices by interpolation on Map array

if keyword_set(NearestNeighbour) then begin
	L = round(L)
	Map = Map( reform(L[0,*]), reform(L[1,*]))
endif else begin
	Map = reform( Interpolate(Map,L[0,*],L[1,*]) )
	Map = (round(Map) > MapMin) < MapMax
endelse							; Map now contains one element for each vertex

MapMin = min(Map,max=MapMax)

; At this point Map scales from Colors[0] to Colors[1]
; Keyword Scale limits then number of colors to only Scale[0] evenly spaced
; in the range Colors[0] to Colors[1]

if n_elements(Scale) ne 0 then begin
	S	= Scale[0] < (MapMax-MapMin+1)
	L   = MapMin+(MapMax-MapMin)/float(S)*(1+findgen(S-1))
	Map = MapMin+(MapMax-MapMin)*(Scale2Levels(Map, L)/float(S-1))
endif

if n_elements(Polygon) ne 0 then begin
	nP = n_elements(Polygon)
	P = 0L		; Initialize P (this element is dropped again later)

	i = -1L
	repeat begin
		i = i+1
		P = [P,Polygon[i+1]]	; Pick up first vertex for each polygon
		i = i+Polygon[i]
		print,i,np
	endrep until i ge nP-1
	P = P[1:*]
	Map = Map[P]				; Extract list of first vertices from full list
								; Map now contains one element for each polygon
endif

return, Map  &  end