| Re: bounding box subsetting for lat/lon [message #69849] |
Tue, 16 February 2010 07:18 |
penteado
Messages: 866
Registered: February 2018
|
Senior Member
Administrator |
|
|
On Feb 16, 4:27 am, mankoff <mank...@gmail.com> wrote:
> But at a minimum I want to specify 4 lat/lon elements that do not
> defined the edges of the box but the four corners, so the roi could be
> a diamond shape (when projected on, say, the /cylinder projection).
> And ideally I could define a path of n points defining a closed region
> and get the indices back that lie within that region.
>
> Anyone know if this exists? I don't see it under the ROI section, and
> searching for bounding box brings up too much noise.
For the generic problem of finding the points inside a region in a
plane, I would use an ROI:
oroi=obj_new('idlanroi',xverts,yverts)
interior=oroi->containspoints(datax,datay)
w=where(interior ne 0,nw)
if (nw gt 0) then begin
datax=datax[w]
datay=datay[w]
dataval=dataval[w]
endif else (...)
obj_destroy,oroi
The trouble with the coordinates being lat/lon is that they are
spherical coordinates. That means if you used that for lon as x and
lat as y, you would actually be dealing with a polygon in a Cartesian
space (cylindrical projection). Which may be what you want, but is not
the same as the closed region in the surface of the sphere that has
these vertices (which would presumably be connected by great circles,
not straight lines). Also, any closed shape on the surface of a sphere
defines two closed regions, depending on which side of the line you
consider the inside.
|
|
|
|
| Re: bounding box subsetting for lat/lon [message #69853 is a reply to message #69849] |
Tue, 16 February 2010 05:04  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
mankoff writes:
> I'm wondering if a routine already exists that can subset a data set
> based on non-square lat/lon limits. I have large 1D arrays of lat,
> lon, and data, and want to do something like this:
>
> roi = WHERE( lat LE a AND lat GE b AND lon LE c AND lon GE d )
> lat = lat[roi] & lon=lon[roi] & data = data[roi]
>
> But at a minimum I want to specify 4 lat/lon elements that do not
> defined the edges of the box but the four corners, so the roi could be
> a diamond shape (when projected on, say, the /cylinder projection).
> And ideally I could define a path of n points defining a closed region
> and get the indices back that lie within that region.
>
> Anyone know if this exists? I don't see it under the ROI section, and
> searching for bounding box brings up too much noise.
The biggest mistake people make in working with map projections
is working in lat/lon space. This is not a space where you
can do much of anything worthwhile in terms of selecting and
manipulating pixels. Life will be MUCH simpler if you work
in Cartesian XY space, where it is quite easy to define your
"box" by either its top, bottom and sides, or by any two opposite
corners. You can use the Map_Proj_Inverse and Map_Proj_Forward
routines to go back and forth between geographic and Cartesian
coordinates.
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
