; Draw the BAS Magic DEM to the screen
set_plot,'z'
res=2000
device,set_res=res*[1,1]
@set_map

; Read the screen into z
z=tvrd()               

; Read in the proposed mask, for the grid
a=get_field('higem-proposed-sh-mask.pp')
la=yc(a,/arr)
lo=xc(a,/arr)

; Convert our lat-lon to device 
xy=convert_coord(lo,la,/to_dev,/data) 
i=where(xy (0,*) ge 0 and xy(0,*) le res and xy(1,*) ge 0 and xy(1,*) le res)

; Now fill in a copy of the original mask with 0's and 1's
; Nb: in z,
;     background (sea) = 255's.
;     plotted (land)   = 0's.
; (depends on colour table?)
b=float(z(xy(0,i),xy(1,i))/255)                                                         
a1=a
a1.bmdi=1.0
a1.data=0.0
a1.data(i)=reform(float(b))                                                  

; Result is now in "a1"

set_plot,'x'