; returns the file X and Y coordinates from the original GPS coordinates

; eg.
;   dem = read_tiff('e:\arolla\arolla.tif',geotiff=geo)
;   xyfile=map2fileXY(dem, geo, 606090,90337)
;Eg.
;   X=[606583.,607037.,607781.]
;   Y=[91342.0,90566.4,90337.4]
;   xyfile=map2fileXY(dem, geo, X,Y, /D)
;
; settin /D display the coordinates on the DEM
;Javier G. Corripio July, 2004

function map2fileXY, dem, geo, xcoord, ycoord, D=D
On_error,2
IF (n_params() LT 4) THEN BEGIN
    message, 'USAGE: Result = map2fileXY(dem, geotiff, xcoord, ycoord, [/D])',/INFO
        RETURN,!values.F_NaN
ENDIF

;;** check if DEM exists and load it  (may read directly....
;    checktiff=FINDFILE(filename)
;    IF (filename NE checktiff[0]) THEN BEGIN
;    filename=DIALOG_PICKFILE(TITLE='Please select DEM file',/MUST_EXIST)
;    WIDGET_CONTROL,demname,SET_VALUE=filename
;    ENDIF
;    z=read_tiff(filename,GEOTIFF=geo)

 ;** DEM information
    upperleftx=geo.modeltiepointtag[3]
    upperlefty=geo.modeltiepointtag[4]
    dl=geo.modelpixelscaletag[0] ;assumes regular grid
    cr=size(dem)
    cols=FIX(cr[1])
    rows=FIX(cr[2])
    Ymin=upperlefty-(dl*rows)
    Xmax=upperleftX+(dl*cols)

    filex = ROUND((xcoord-upperleftx)/dl)
    filey = ROUND((upperlefty-ycoord)/dl)

    IF (max(filex) gt cols) or (max(filey) Gt rows) or $
       (min(filex) lt 0) or (min(filey) LT 0) THEN BEGIN
       message, 'Coordinates are outside DEM range',/INFO
        RETURN,!values.F_NaN
    ENDIF
    dimxy = size(filex,/dimensions)
    dimxy=dimxy[0]
    filexy = fltarr(2,dimxy)
    filexy[0,*]=filex
    filexy[1,*]=filey

    IF KEYWORD_SET(D) THEN BEGIN
       demimg=dem
       for i=0, dimxy-1 do begin
           demimg[filex[i]-3:filex[i]+3,filey[i]]=0
           demimg[filex[i],filey[i]-3:filey[i]+3]=0
        endfor
       display,demimg
    ENDIF



RETURN,filexy


end