Re: weird behavior of Triangulate [message #81355 is a reply to message #81286] |
Tue, 11 September 2012 12:53   |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Yngvar Larsen writes:
> ;; Create IDL coordinate mapper for polar stereographic grid.
> ;; See http://nsidc.org/data/polar_stereo/ps_grids.html
> re = 6378273d0
> e = 0.081816153d0
> rp = sqrt((1-e^2)*re^2)
> map = map_proj_init('polar stereographic', /gctp, $
> center_latitude=70, $
> center_longitude=315, $
> semimajor_axis=re, $
> semiminor_axis=rp)
>
> ;; Define output grid.
> dx = 25d3
> dy = 25d3
> xmin = -3.85d6
> ymin = -5.35d6
> nxout = 304L
> nyout = 448L
> mapx = xmin + dx*dindgen(nxout)
> mapy = ymin + dy*dindgen(nyout)
>
> mapx = mapx[*,lindgen(nyout)]
> mapy = transpose(mapy[*,lindgen(nxout)])
>
> ;; Calculate output map grid values in input coordinates
> latlon = map_proj_inverse(mapx, mapy, map_structure=map)
I don't understand this step. In your explanation of the
"right way" you say:
"Start with a regular grid in the _output_ domain.
Calculate where these grid points are located in
the input domain."
I would have thought this means "use the map structure
of the input data" in this case. And, yet, you are using
the map structure of the output data grid. Yes, the
points will all fall within the input domain, since
that is a global domain. But, what about when the
input domain is NOT global, so that some of the points
are inside the output domain, but some are outside, too?
I can see that your program works, and it is EXTREMELY
fast. I just can't see *why* it works yet. :-)
Thanks for your help. I'm learning a lot! :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|