Matching Lats and Lons from two arrays

Paolo,

thanks you were spot on. That bought me another x10 speed improvemnt
in that routine.

Thanks,

Brian

--------------------------------------------------------------------------
Brian Larsen
Boston University
Center for Space Physics
http://people.bu.edu/balarsen/Home/IDL

Forgive my trespass but why is nobody suggesting the use of the
DISTANCE_MEASURE or MAP_2POINTS functions of IDL? Or is the problem
only with the length of the loop?

Or did I not get the problem right?

Gaurav

On Aug 26, 11:47*am, wilsona <awils...@bigred.unl.edu> wrote:
> I have 2 seperate arrays of Latittudes and Longitudes.
> CS_LATLON(0,4607) is one latitude array and dlat(192,139) is the
> other
> CS_LATLON(1,4607) is one longitude array and dlon(192,139) is the
> other.
> I want to index through each element in both CS_LATLON arrays and
> find
> which point(s) in the dlat and dlong arrays are closest.
> I tried using nested loops but that gave me 12 million+ loops which
> is
> too many for my liking. *I now am trying search2d
> NUM_PNTS *= N_ELEMENTS(CS_LATLON(0, *)) - 1
>
> * FOR J = 0, NUM_PNTS DO BEGIN
> * * * CLOSE_LATS = SEARCH2D(dlat, 0, 0, CS_LATLON(0,J),
> CS_LATLON(0,J), INCREASE=0.5, * *$
> * * * * * * * * * * * * * * DECREASE=0.5)
> * * * lat1 * = CS_LATLON(0,J) * PI / 180.0
> * * * FOR K = 0, CLOSE_LATS DO BEGIN
> * * * * * * * lat2 * = dlat(K) * PI / 180.0
> * * * * * * * d_long = CS_LATLON(1,J) - dlon(K)) * PI / 180.0
> * * * * * * * DISTANCE * = 10800.0 / PI * acos(sin(lat1) * sin(lat2)
> +
> cos(lat1) * * * * $
> * * * * * * * * * * * * * * * * * * *cos(lat2) * cos(d_long))
> * * * *ENDFOR ; K
> * ENDFOR ; J
> This is not working they way I would like. *Any suggestions on this
> would be greatly appreciated.

You might find WITHINSPHRAD in JBIU useful:

http://astroconst.org/jbiu/jbiu-doc/...hinsphrad.html

-Jeremy.

I should point out that the MATCH_2D routine mentioned in the above
thread and at:

http://www.dfanning.com/code_tips/matchlists.html

assumes euclidean distance measures apply, which of course is *not*
strictly correct on a sphere. For points of latitude and longitude
(or ra/dec), this will work for small patches of earth (or sky), but
if you are near either pole, or cover appreciable areas, the 2D
distance measure will falter.

You can "fake" a "pretty good" euclidean distance measure by pre-
multiplying all longitudes by the cos of latitude. This works if the
latitude range is not overly large, and if you don't cross any poles.
Do do this correctly would require replacing the d=dx^2 + dy^2 with:

d = 2 asin( sqrt( sin(ddec/2)^2 + cos(dec1)cos(dec2)sin(dra/2)^2 ) )

which is obviously much costlier to evaluate. More importantly,
you'll need to think a bit about whether the stock "platte carre"
projection (i.e. (ra, dec) -> (x,y)) is wasteful or not for your
distribution of points on the sphere, or adversely affects your
matching radius criterion. Adopting another projection prior to
binning and running the histogram might offer gains in efficiency and
more "uniform" distance performance (the issue being that square bins
in a projected sphere do not have constant area). It should be
straightforward to "wrap-around" the poles.

JD