; 
; grid the data and display an image
; then you can save as .png from iimage
;SYNTAX:istat=griddata(xyzIn,radErr,img,xout,yout,tr=tr,useTr=useTr,xyzOut=xyzOut,minpnts=minpnts,$
;                       searchellipse=searchellipse,gridLen=gridlen,axisF=axisF,$
;                       rotateGrid=rotateGrid,verb=verb)
;-------------------------------------------------------
;ARGS:
; xyzIn[3,n]: dbl  xyz points to grid (in meters)
; radErr[n]:dbl  radial error for each xyzM point. can be meters, cm, whatever you want
;
;	`
;KEYWORDS:
;minPnts    : long # of points to use around grid point. def=40
;searchEllipse: meters .. limit search for points to searchEllipse. def is .5 meter
;                         if no data within this region then mark as no data.
;                         The search ellipse is always specified in meters.
;gridlen    : long    number of points in the gridded data. def=2000 
;axisF      :     if set then convert the x,y values from meters to feet.
;                 The radialErr is not modified.
;rotateGrid : dble  angle(deg) to rotate the grid before triangulation/gridding.
;                   (note: if rotateGrid is not 0 and /usetr is supplied with tr then
;                   the tr must have been computed with the same rotateGrid value)
;                   The rotation is CW about the z axis.
;useTr      :     if set then use tr=tr as the triangulation.Don't bother to compute it
;                 If not set then compute tr and  optionally return it if tr= keyword is present
;tr         :{}   triangulation of xyzM
;verb       :     if set then print progress of the processing
;RETURNS:
;istat   : int       0 - error, not generated
;                    1 - data gridded ok
;img[m,m]: dblarr    gridded data
;xyzOut[3,n]:dblarr   if rotateGrid or axisF set then this will be the xyz grid after the
;                    changes.
;xout[m] : dblarr    xvalues of the grid
;yout[m] : dblarr    yvalues of the grid
;tr               triangulation (can be input an output)
;
;DESCRIPTION:
; 	Take the xyzM locations and the radErrM arrays and generate a gridded image.
;The data in radErrM should already have been limited to the values to used
; eg (say +/- 5cm error). 
;  if rotateGrid is supplied then rotate about about the z axis (warning: if /usetr is set,
;then the supplied tr must have been computed with the same rotation angle.
;
;  triangulate the data (if /usetr is not set) and then grid the data
;using the inverse distance grid.
; the processing order is:
;  1. if axisF make a copy of xyzM[] -> xyz[] and then scale xyz from meters to feet
;  2. if rotateAxis then
;     - if copy  xyzM-> xyz if not already done
;     - rotate xyz by the requested angle
;  3. if usetr is no set, then triangulate  xyz(or xyzM if not yet copied)
;     - if you pass in tr and use it, then it must have been originally computed with the 
;       same rotateAxis and axisF. If not then  tr will not correspond to the modified data.
;  4. grid the data.
;  The game we play with no data
;
; by def we have the gridding routine put -999 in nodata
; after gridding
; 1. we find the (min,max) values (excluding -999)
;    we replace -999 with minvalue - (.001*(max-min)) 
; 2. you can then search for the in value of img and that will be no data
;    and the missingdata will not take up much of the lkup table
;-
function griddatap50,xyzIn,radErr,img,xout,yout,usetr=usetr,tr=tr,xyzOut=xyzOut,minpnts=minpnts,$
					    axisF=axisF,rotateGrid=rotateGrid ,verb=verb,$
						searchellipse=searchellipse,gridLen=gridlen
;
	ftoM=.3048d	
	mTof=1./ftom
	useTr=(keyword_set(usetr))?useTr:0
	if (useTr ne 0) then begin
		if (n_elements(tr) eq 0) then begin
			print,"-->useTr requires you to input tr"
		    return,-1
		endif
	endif
	minPnts=(n_elements(minpnts) eq 1)?minpnts:40l	
	searchEllipseLoc=(n_elements(searchElipse) > 0)?searchEllipse:.5
	gridLen=(n_elements(gridLen) eq 0)?2000:gridLen
	missingData=-999.
	axisF=(n_elements(axisF) eq 0)?0:axisF
	rotateGrid=(n_elements(rotateGrid) eq 0)?0.:rotateGrid
	copyDone=0
	if rotateGrid ne 0. then begin
		xyzOut=rotvec(xyzIn,rotateGrid,axis=3,/usedouble)
		copyDone=1
	endif
	if keyword_set(axisF) then begin
		xyzOut=(copyDone)?xyzOut*mtof:xyzIn*mtof
		copyDone=1
	endif
;
; triangulate
;	
	if not useTr then begin
		; if we rotate to align the cables, triangulate was failing with colinear points
        ; setting the tolerance seems to fit it
		if keyword_set(verb) then print,"Start triangulate:" + systime(0)
	    if copyDone then begin
		    tol=(rotateGrid ne 0.)?max(xyzout)*1d-12:0d
			triangulate,reform(xyzOut[0,*]),reform(xyzOut[1,*]),tr,tol=tol
		endif else begin
			triangulate,reform(xyzIn[0,*])*mtof,reform(xyzIn[1,*])*mtof,tr
		endelse
	endif
;
; grid
;
	scl=(axisF)?mtof:1.
	searchEllipseLoc*=scl
	len=300*scl
	xout=(findgen(gridLen)/(1.*gridLen) - .5)*len
	yout=xout
	if keyword_set(verb) then print,"Start gridding: " + systime(0)
	if copyDone then begin
		img=griddata(reform(xyzOut[0,*]),reform(xyzOut[1,*]),radErr,$
			min_points=minpnts,/inverse_distance,triangles=tr,xout=xout,yout=yout,/grid,missing=missingdata,$
			search_ellipse=searchEllipseLoc)
	endif else begin
		img=griddata(reform(xyzIn[0,*]),reform(xyzIn[1,*]),radErr,$
			min_points=minpnts,/inverse_distance,triangles=tr,xout=xout,yout=yout,/grid,missing=missingdata,$
			search_ellipse=searchEllipse)
	endelse
	
;
; 	switch nodata to minvall - .001(max-min)
;   see if any missing data
;
	jj=where(img eq  missingdata,cntm)
	if cntm gt 0 then begin
;       compute min,max to scale
		ii=where(img ne missingdata,cnt)
		minv=min(img[ii],max=maxv)
		img[jj]=minV - (maxV-minV)*.001
	endif
	if keyword_set(verb) then print,"griddatap50 done: " + systime(0)
	return,1
end