comp.lang.idl-pvwave archive
Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado

Home » Public Forums » archive » Re: Newbie needs help...
Show: Today's Messages :: Show Polls :: Message Navigator
E-mail to friend 
Switch to threaded view of this topic Create a new topic Submit Reply
Re: Newbie needs help... [message #23125] Thu, 11 January 2001 05:29 Go to next message
Bruce Bowler is currently offline  Bruce Bowler
Messages: 128
Registered: September 1998
Senior Member
David Fanning wrote:
>
> Bruce Bowler (bbowler@bigelow.org) writes:
>
>> I'm willing to try this "thinking outside the box" thing for a while,
>> but I can't see the box.
>
> There is no box. (But see more on spiritual development, below.)

Then what did I trip over last night on my way to the
toil<del><del><del><del>reading room?

Fortunately, I think I landed on the outside, and I didn't hurt
myself...

As I hit the ground, I realized that displaying the data and extracting
the data value at some lat/lon are 2 entirely different processes. I
can use Liam's image_map to display it and came up with clever (but as
yet untested) way to extract the data.

Given a target lat/lon and BHAlat and BHAlon, how about (in pseudo-code)

possiblelats = where(BHAlat eq lat{+/- some epsilon})
possiblelons = where(BHAlon eq lon{+/- some epsilon})
possiblevalues = intersection(possiblelats,possiblelons)

if number of possiblevalues is between 1 and 10, printout the data
otherwise adjust epsilon either up or down and try again.

Does that sound like it ought to work (and in some time less than a
glacial epoch) ?

Bruce
Re: Newbie needs help... [message #23133 is a reply to message #23125] Wed, 10 January 2001 14:55 Go to previous messageGo to next message
James Kuyper is currently offline  James Kuyper
Messages: 425
Registered: March 2000
Senior Member
"Liam E. Gumley" wrote:
...
> You mentioned that you are using MODIS data at 1 km resolution (1354 x
> 2030 pixels pre granule). One approach we use is to define a global
> equal area grid; we happen to use 25 km x 25 km grid cells. We loop over
> each 1 km pixel, and based on it's lat/lon, we accumulate the following
> statistics in each grid cell:
>
> Number of observations,
> Sum of observations,
> Minimum observation,
> Maximum observation.
>
> From these, we can compute mean and standard deviation. it is quite
> straightforward to then resample the equal area grid to an equal angle
> grid which can be visualized in IDL.

One thing to keep in mind about MODIS data is that you don't need to
grid it yourself. For most data sets, if MODIS has a pixel-oriented
product, it also has one or more corresponding products which are
gridded by latitude and longitude. For instance, the MOD10_L2 product
reports on snow coverage for each pixel, with 500m resolution at nadir;
each file covering 5 minutes of data. MOD10L2G, on the other hand,
summarizes an entire day's worth of data gridded according to one of
several possible different map projections, with each file covering a
tile that is nominally 10 deg by 10 degrees, with a nominal grid
resolution of 500m.
Re: Newbie needs help... [message #23134 is a reply to message #23133] Wed, 10 January 2001 14:46 Go to previous messageGo to next message
Craig Markwardt is currently offline  Craig Markwardt
Messages: 1869
Registered: November 1996
Senior Member
"Liam E. Gumley" <Liam.Gumley@ssec.wisc.edu> writes:
> You mentioned that you are using MODIS data at 1 km resolution (1354 x
> 2030 pixels pre granule). One approach we use is to define a global
> equal area grid; we happen to use 25 km x 25 km grid cells. We loop over
> each 1 km pixel, and based on it's lat/lon, we accumulate the following
> statistics in each grid cell:
>
> Number of observations,
> Sum of observations,
> Minimum observation,
> Maximum observation.
>
> From these, we can compute mean and standard deviation. it is quite
> straightforward to then resample the equal area grid to an equal angle
> grid which can be visualized in IDL.

Generally speaking I am a scientist who looks *away* from the earth
not towards it, but with satellite-based astronomy there are still
times where one has to worry about where the satellite is.

I am occassionally charged with constructing various instrument
housekeeping data. Typically we receive multiple measurements at the
same geographic point and I needed to average them. I wrote this
little routine which constructed the mean and variance of the data.
It's called GEOSPLAT, but there's really nothing "GEO" about it. It
simply constructs a histogram of LAT/LON data points, weighted by the
measurement values. Thankfully our satellite is in a low inclination
orbit, so I never really deal with the poles. This routine just makes
a 2-d histogram regularly sampled in LON/LAT space.

This routine requires a specially modified version of HIST_2D called
HIST_2DR, but don't worry, all I did was add the REVERSE_INDICES
keyword to HIST_2D. I can forward that if needed.

This application of REVERSE_INDICES is actually not the most
efficient, since it loops through every cell in the output array.
There is another way to do it with far fewer iterations. If you have
very large image, then one should shift to that.

Craig

function geosplat, lon, lat, val, expo, xrange=xrange0, yrange=yrange0, $
xbinsize=xbinsize, ybinsize=ybinsize, minsamp=minsamp, $
variance=hh2, exposure=ee

forward_function arg_present, hist_2dr

if n_elements(xrange0) LT 2 then xrange0 = [0,359.9999]
if n_elements(yrange0) LT 2 then yrange0 = [-90,89.9999]
xrange = xrange0(0:1)
yrange = yrange0(0:1)
if n_elements(xbinsize) LT 1 then xbinsize = 1
if n_elements(ybinsize) LT 1 then ybinsize = 1
if n_elements(minsamp) LT 1 then minsamp = 1
if n_elements(expo) LT n_elements(val) then expo = val*0 + 1.
if arg_present(hh2) then dohh2 = 1 else dohh2 = 0

hh = hist_2dr(lon, lat, reverse=rr0, $
min1=xrange(0), max1=xrange(1), bin1=xbinsize(0), $
min2=yrange(0), max2=yrange(1), bin2=ybinsize(0))
hh = hh * (val(0)*0)
ee = hh
if dohh2 then hh2 = hh
for i = 0, n_elements(hh)-1 do begin
if rr0(i+1) - rr0(i) GT minsamp(0) then begin
ex = expo(rr0(rr0(i):rr0(i+1)-1))
vv = val(rr0(rr0(i):rr0(i+1)-1))
hh(i) = total(vv*ex)
ee(i) = total(ex)
if dohh2 then hh2(i) = total(vv^2*ex)
endif
endfor

hh = hh / (ee>1e-20)
hh2 = hh2 / (ee>1e-20)
hh2 = hh2-hh^2

return, hh
end



--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
Re: Newbie needs help... [message #23135 is a reply to message #23134] Wed, 10 January 2001 12:49 Go to previous messageGo to next message
Liam E. Gumley is currently offline  Liam E. Gumley
Messages: 378
Registered: January 2000
Senior Member
David Fanning wrote:
> (Where *is* Liam, anyway!?)

I'm right here, trying to get my semi-annual research report written...

> I guess if I were in your shoes, I might try to
> reformat the part of the image I was interested
> in onto a regular lat/lon grid. I'd probably try
> to use TRIANGULATE and TRIGRID first, to see if
> that would work.

In my experience, TRIANGULATE and TRIGRID work just fine when the input
arrays are small (say a few hundred rows and columns). But when you have
satellite images with thousands of rows and columns, TRIANGULATE and
TRIGRID just don't work.

> There are probably other ways to warp the image
> onto the map, but there is no getting around the
> need to do it. What you would prefer, I'm sure, is
> to have a map warped to the image. But this is NOT
> possible in IDL.

Displaying a satellite image on a map projection is one problem, and my
imagemap procedure provides one solution.

However resampling a satellite image to a regular lat/lon grid in such a
way that data values are preserved is a different problem, and there is
no easy solution in IDL. Typically, the approach is to define your
lat/lon grid (in some map projection) and then for each map grid cell,
find the closest matching location in the satellite image (nearest
neighbor). Refinements include bilinear or cubic spline interpolation.
However this algorithm is not easy to implement efficiently in IDL,
because of the penaly associated with loops. Most people implement a
resampling algorithm in some other language (e.g. C or FORTRAN), and
then read the results in IDL.

You mentioned that you are using MODIS data at 1 km resolution (1354 x
2030 pixels pre granule). One approach we use is to define a global
equal area grid; we happen to use 25 km x 25 km grid cells. We loop over
each 1 km pixel, and based on it's lat/lon, we accumulate the following
statistics in each grid cell:

Number of observations,
Sum of observations,
Minimum observation,
Maximum observation.

From these, we can compute mean and standard deviation. it is quite
straightforward to then resample the equal area grid to an equal angle
grid which can be visualized in IDL.

If you are interested in more details, contact me directly.

Now back to my semi-annual report (sigh).

Cheers,
Liam.
http://cimss.ssec.wisc.edu/~gumley

PS: Here's a couple of MODIS images created in IDL:
http://earthobservatory.nasa.gov/Newsroom/NewImages/images.p hp3?img_id=4549
http://earthobservatory.nasa.gov/Newsroom/NewImages/images.p hp3?img_id=4383
Re: Newbie needs help... [message #23136 is a reply to message #23135] Wed, 10 January 2001 12:18 Go to previous messageGo to next message
davidf is currently offline  davidf
Messages: 2866
Registered: September 1996
Senior Member
Bruce Bowler (bbowler@bigelow.org) writes:

> karmic confusion reigns supreme...

Now we *are* getting somewhere... :-)

> I have this big huge array (here
> after referred to as a BHA) of data (BHAd) an equally sized BHA of
> latitudes (BHAlat) of each point in BHAd and another BHA of longitudes
> (BHAlon) of each point in BHAd. NB, row x of BHAlat is not constant,
> neither is column y of BHAlon
>
> I understand map_set, I understand map_grid, I understand map_continents
> (well, not in the tao-ist sense of the word). map_image is where I'm
> getting confused... how does map_image know about BHAlat and BHAlon
> when they are not inputs to the process?

Ah, this is what IDL doesn't do. Liam does this
by taking each one of your pixels (I'm paraphrasing
here) and putting it over on your map one by one.
Then he sees that he has a very holey map (pun
intended in this case), so he "fills it in" with
some kind of dilation thing, if I remember correctly.
(Where *is* Liam, anyway!?) That's why he doesn't
want you to trust the numbers. He has had to fudge
them a bit. But they look great!

I guess if I were in your shoes, I might try to
reformat the part of the image I was interested
in onto a regular lat/lon grid. I'd probably try
to use TRIANGULATE and TRIGRID first, to see if
that would work. You will jimmy the numbers, of
course, but no way around that if you want to use
IDL. Then, I would put the regular lat/lon image
on the map using Map_Image.

There are probably other ways to warp the image
onto the map, but there is no getting around the
need to do it. What you would prefer, I'm sure, is
to have a map warped to the image. But this is NOT
possible in IDL.

Then, I would get the specified value from the
regular lat/lon grid, if I couldn't figure out
a way to get a better number from my original
array. Remember that the display *IS* only a
display. You are not really looking at your data.
All you want to do is have the display *represent*
your data in a more or less faithful way. I think
Liam's solution is as good as you are likely to
get for actually displaying your data.

Once you have the display, getting the actual value
out of your data set shouldn't be terribly hard.
You might need some kind of look-up table for
indexing, but this should be fairly straightforward.

Cheers,

David

P.S. Let's just say I really wish Liam would step
in here. :-(

--
David Fanning, Ph.D.
Fanning Software Consulting
Phone: 970-221-0438 E-Mail: davidf@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Toll-Free IDL Book Orders: 1-888-461-0155
Re: Newbie needs help... [message #23137 is a reply to message #23136] Wed, 10 January 2001 11:44 Go to previous messageGo to next message
Bruce Bowler is currently offline  Bruce Bowler
Messages: 128
Registered: September 1998
Senior Member
karmic confusion reigns supreme...

David Fanning wrote:
>
> Bruce Bowler (bbowler@bigelow.org) writes:
>
>> I'm willing to try this "thinking outside the box" thing for a while,
>> but I can't see the box.
>
> There is no box. (But see more on spiritual development, below.)
>
>> Now that I have my image mapped on to a lat/lon grid (see, I'm learning
>> already :-), how to I access the data by lat/lon?
>
> You are going to love this! :-)
>
> You have set up the map coordinate space with Map_Set.
> You have placed your image on the map with Map_Image
> (or something similar). You put your map grid and
> continental boundary on your map with Map_Grid
> and Map_Continents.

now we maybe be getting somewhere... I have this big huge array (here
after referred to as a BHA) of data (BHAd) an equally sized BHA of
latitudes (BHAlat) of each point in BHAd and another BHA of longitudes
(BHAlon) of each point in BHAd. NB, row x of BHAlat is not constant,
neither is column y of BHAlon

I understand map_set, I understand map_grid, I understand map_continents
(well, not in the tao-ist sense of the word). map_image is where I'm
getting confused... how does map_image know about BHAlat and BHAlon
when they are not inputs to the process?

To add a little to the confusion, and to give some scope to the problem,
BHAd is an array [1354,2030] and covers lat/lon
[34.9,-78.33,56.58,-41.68]. I'm only interested in lat/lon
[41,-71,45,-66] (it's a MODIS swath but the data's an experimental
product so none of my "normal" MODIS tools work), we're doing sea-truth
work for this product in the Gulf of Maine.

> Now, you get lat/lon value from the user by just
> having them click on the map! Too easy!
>
> Cursor, lon, lat, /Data

no, they're going to tell me the lat and lon from a data file, but that
part I can handle.

> Now, what you do next depends on you. If you
> have an image data set in which each pixel
> has an associated lat/lon coordinate, you can
> go pull out the closest pixel value from that
> data set.

I suspect once I figure out the map_image part, this part MIGHT fall out
on it's own...

> If you don't have such a data set, you might
> have to get the value from the warped image.
> That value, of course, was created by smushing
> (technical term) several real values together
> in the warping process. It may not be what you
> want. (How this is done is fairly complicated.
> I'd explain it, but I'm pretty sure it's not
> what you want to do anyway.)

Ah shucks, go ahead, explain it. It can't make my head hurt anymore
than it already does...
Re: Newbie needs help... [message #23139 is a reply to message #23137] Wed, 10 January 2001 09:25 Go to previous messageGo to next message
davidf is currently offline  davidf
Messages: 2866
Registered: September 1996
Senior Member
Bruce Bowler (bbowler@bigelow.org) writes:

> I'm willing to try this "thinking outside the box" thing for a while,
> but I can't see the box.

There is no box. (But see more on spiritual development, below.)

> Now that I have my image mapped on to a lat/lon grid (see, I'm learning
> already :-), how to I access the data by lat/lon?

You are going to love this! :-)

You have set up the map coordinate space with Map_Set.
You have placed your image on the map with Map_Image
(or something similar). You put your map grid and
continental boundary on your map with Map_Grid
and Map_Continents.

Now, you get lat/lon value from the user by just
having them click on the map! Too easy!

Cursor, lon, lat, /Data

Now, what you do next depends on you. If you
have an image data set in which each pixel
has an associated lat/lon coordinate, you can
go pull out the closest pixel value from that
data set.

If you don't have such a data set, you might
have to get the value from the warped image.
That value, of course, was created by smushing
(technical term) several real values together
in the warping process. It may not be what you
want. (How this is done is fairly complicated.
I'd explain it, but I'm pretty sure it's not
what you want to do anyway.)

>> if you persist I can recommend "getting used
>> to disappointment" as the appropriate spiritual practice.
>>
>> Cheers,
>
> Those two sentiments don't seem to go together :-(

Oh, I don't know. I find it an easier spiritual
practice than "do the dishes, again", which is what
I usually do in my search for enlightenment. :-)

Cheers,

David
--
David Fanning, Ph.D.
Fanning Software Consulting
Phone: 970-221-0438 E-Mail: davidf@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Toll-Free IDL Book Orders: 1-888-461-0155
Re: Newbie needs help... [message #23140 is a reply to message #23139] Wed, 10 January 2001 08:33 Go to previous messageGo to next message
Bruce Bowler is currently offline  Bruce Bowler
Messages: 128
Registered: September 1998
Senior Member
David Fanning wrote:
>
> Bruce Bowler (bbowler@bigelow.org) writes:
>
>> I've lurked for a while and now it's time to ask for help...
>
> Uh, oh. Trouble here. Newbies and de-cloaked lurkers
> always ask the most difficult questions. :-(

It get's worse :-) I've been programming in THE classical programming
language (FORTRAN roolz!) for more the 2/3's of my life (I'm 45 now, you
do the math :-). The first machine I ever wrote a program for was the
size of a small bus and had a whopping 16K (yes, that's a K) of memory.
But let's not start swapping war stories here...

>> I have some satellite data that I need to a) display a portion of on the
>> globe and b) extract the data at various positions.
>
> What you want to do--slap a lat/lon grid on top of
> an image--is reasonable. It's just that IDL is not
> the software you should be using to do it. :-)
>
> IDL works the other way around. It allows you to
> slap an image on top of a lat/lon grid. Sometimes
> the result is similar, but if you hold pixel values
> to be sacrosanct, then you are almost always disappointed
> in what you can do with IDL.

I prefer to hold pixels are arm's length (remember, I'm 45 and I've been
staring at CRT's for a good number of those years, my eyesight ain't
what it used to be :-)

I'm willing to try this "thinking outside the box" thing for a while,
but I can't see the box.

Now that I have my image mapped on to a lat/lon grid (see, I'm learning
already :-), how to I access the data by lat/lon?

> if you persist I can recommend "getting used
> to disappointment" as the appropriate spiritual practice.
>
> Cheers,

Those two sentiments don't seem to go together :-(
Re: Newbie needs help... [message #23142 is a reply to message #23140] Wed, 10 January 2001 08:08 Go to previous messageGo to next message
davidf is currently offline  davidf
Messages: 2866
Registered: September 1996
Senior Member
Bruce Bowler (bbowler@bigelow.org) writes:

> I've lurked for a while and now it's time to ask for help...

Uh, oh. Trouble here. Newbies and de-cloaked lurkers
always ask the most difficult questions. :-(

> I have some satellite data that I need to a) display a portion of on the
> globe and b) extract the data at various positions.
>
> Towards that end, I've snarfed up image_map.pro from Liam Gumley's site
> and regrid.pro and tvim.pro from ESRG. They both do sort of what I want
> but neither gets me "all the way there".
>
> Liam's does a nice job of displaying the image and letting me plop a
> grid and continent map on top but he's got this big fat caveat in the
> code that says "This procedure was designed for display purposes *only*"
> (emphasis his) and some bits about dire consequences if used otherwise.
>
> regrid does a good job of "regularizing" the data and tvim makes a nice
> plot with scales and color bars and all that nice stuff (without the
> "dire consequences" :-) but when I try and slap a grid or continental
> blobs on the image, they puke complaining that "the current device must
> have mapping coordinates"...
>
> In either case, I'm left with the problem of mapping a lat/lon that the
> user inputs to a pixel on the screen (or an element in an array, as the
> case may be).

What you want to do--slap a lat/lon grid on top of
an image--is reasonable. It's just that IDL is not
the software you should be using to do it. :-)

IDL works the other way around. It allows you to
slap an image on top of a lat/lon grid. Sometimes
the result is similar, but if you hold pixel values
to be sacrosanct, then you are almost always disappointed
in what you can do with IDL.

Liam offers his big fat caveat with good reason. He figured out
a clever way of making the image display *look* right, but
he is rightfully offering no guarantees about the numbers.

Liam can fill you in on the details of this better than
I can, but if you persist I can recommend "getting used
to disappointment" as the appropriate spiritual practice.

Cheers,

David
--
David Fanning, Ph.D.
Fanning Software Consulting
Phone: 970-221-0438 E-Mail: davidf@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Toll-Free IDL Book Orders: 1-888-461-0155
Re: Newbie needs help... [message #23221 is a reply to message #23125] Thu, 11 January 2001 07:53 Go to previous message
Liam E. Gumley is currently offline  Liam E. Gumley
Messages: 378
Registered: January 2000
Senior Member
Bruce Bowler wrote:
> As I hit the ground, I realized that displaying the data and extracting
> the data value at some lat/lon are 2 entirely different processes. I
> can use Liam's image_map to display it and came up with clever (but as
> yet untested) way to extract the data.
>
> Given a target lat/lon and BHAlat and BHAlon, how about (in pseudo-code)
>
> possiblelats = where(BHAlat eq lat{+/- some epsilon})
> possiblelons = where(BHAlon eq lon{+/- some epsilon})
> possiblevalues = intersection(possiblelats,possiblelons)
>
> if number of possiblevalues is between 1 and 10, printout the data
> otherwise adjust epsilon either up or down and try again.

I've written something similar for the MODIS Airborne Simulator (MAS),
which you could probably adapt for MODIS by tuning the epsilon values.
It works reasonably efficiently for small numbers of pixels. A couple of
other routines are required:

setintersection.pro from RSI:
http://www.dfanning.com/tips/set_operations.html

compass.pro from ESRG:
http://www.astro.washington.edu/deutsch-bin/idllibsrch?keywo rd=compass

This routine is appropriate for finding a few lat/lon locations at a
time. However it would not be very effective for overlaying coastline
lat/lon vectors on an image.

Cheers,
Liam.
http://cimss.ssec.wisc.edu/~gumley

;-----------------------------------------
PRO MAS_LOCATE, SLAT, SLON, LAT, LON, X, Y

;+
; PURPOSE:
; Locate a given lat/lon in a MAS image.
;
; INPUT:
; SLAT Latitude to locate (deg)
; SLON Longitude to locate (deg)
; LAT Array of MAS latitude values (deg)
; LON Array of MAS longitude values (deg)
;
; OUTPUT:
; X Pixel number closest to the given lat/lon (-1 if not found)
; Y Line number closest to the given lat/lon (-1 if not found)
;
; REVISED:
; Liam.Gumley@ssec.wisc.edu
; $Id: mas_locate.pro,v 1.2 1999/10/29 16:26:49 gumley Exp $
;-

;- Check arguments

if n_params() ne 6 then message, 'Usage: MAS_LOCATE, SLAT, SLON, LAT,
LON, X, Y'
if n_elements(slat) eq 0 then message, 'SLAT is undefined'
if n_elements(slon) eq 0 then message, 'SLON is undefined'
if n_elements(lat) eq 0 then message, 'LAT is undefined'
if n_elements(lat) eq 0 then message, 'LON is undefined'
if size(lat, /n_dim) ne 2 then message, 'LAT is not a 2D array'
if size(lon, /n_dim) ne 2 then message, 'LON is not a 2D array'
if arg_present(x) eq 0 then message, 'X cannot be modified'
if arg_present(y) eq 0 then message, 'Y cannot be modified'

;- Set default return values

x = -1L
y = -1L

;- Check that lat/lon lies within the array min/max

latmin = min(lat, max=latmax)
lonmin = min(lon, max=lonmax)
if (slat lt latmin) or (slat gt latmax) or $
(slon lt lonmin) or (slon gt lonmax) then return

;- Find array elements close to the lat/lon

latindex = where(abs(lat - slat) lt 0.001, latcount)
lonindex = where(abs(lon - slon) lt 0.001, loncount)
if (latcount lt 1) or (loncount lt 1) then return

;- Find the intersecting elements of the arrays

result = setintersection(latindex, lonindex)
if (result[0] eq -1) then return

;- Compute the distance from the lat/lon to the array elements

compass, slat, slon, lat[result], lon[result], range, azimuth

;- Find the array element closest to the lat/lon

minrange = min(range, minindex)

;- Convert the 1D array element index to x/y

dims = size(lat, /dim)
x = result[minindex] mod dims[0]
y = result[minindex] / dims[0]

END
;-----------------------------------------
Re: Newbie needs help... [message #23223 is a reply to message #23125] Thu, 11 January 2001 06:45 Go to previous message
Paul van Delst is currently offline  Paul van Delst
Messages: 364
Registered: March 1997
Senior Member
Bruce Bowler wrote:
>
> David Fanning wrote:
>>
>> Bruce Bowler (bbowler@bigelow.org) writes:
>>
>>> I'm willing to try this "thinking outside the box" thing for a while,
>>> but I can't see the box.
>>
>> There is no box. (But see more on spiritual development, below.)
>
> Then what did I trip over last night on my way to the
> toil<del><del><del><del>reading room?
>
> Fortunately, I think I landed on the outside, and I didn't hurt
> myself...
>
> As I hit the ground, I realized that displaying the data and extracting
> the data value at some lat/lon are 2 entirely different processes. I
> can use Liam's image_map to display it and came up with clever (but as
> yet untested) way to extract the data.
>
> Given a target lat/lon and BHAlat and BHAlon, how about (in pseudo-code)
>
> possiblelats = where(BHAlat eq lat{+/- some epsilon})
> possiblelons = where(BHAlon eq lon{+/- some epsilon})
> possiblevalues = intersection(possiblelats,possiblelons)
>
> if number of possiblevalues is between 1 and 10, printout the data
> otherwise adjust epsilon either up or down and try again.
>
> Does that sound like it ought to work (and in some time less than a
> glacial epoch) ?
>
> Bruce

It seems to me that the problem is one where the transformation from array lat/lon to map lat/lon
(for viewing as a contiguous eye-candy type of image) is linear, but the reverse is decidedly
non-linear. Your nearest neighbour search above seems like a workable solution since your
referencing the original lat/lon array.

paulv

--
Paul van Delst A little learning is a dangerous thing;
CIMSS @ NOAA/NCEP Drink deep, or taste not the Pierian spring;
Ph: (301) 763-8000 x7274 There shallow draughts intoxicate the brain,
Fax: (301) 763-8545 And drinking largely sobers us again.
Email: pvandelst@ncep.noaa.gov Alexander Pope.
Re: Newbie needs help... [message #23224 is a reply to message #23125] Thu, 11 January 2001 06:59 Go to previous message
davidf is currently offline  davidf
Messages: 2866
Registered: September 1996
Senior Member
Bruce Bowler (bbowler@bigelow.org) writes:

> As I hit the ground, I realized that displaying the data and extracting
> the data value at some lat/lon are 2 entirely different processes. I
> can use Liam's image_map to display it and came up with clever (but as
> yet untested) way to extract the data.

*Now* you are on the path! :-)

Too many people, I think, fail to make the distinction
between the *display* of their data, and the data itself.
They are almost never the same thing. If you can keep
the distinction clear, and always have a way to get
to the real data in your program, then you have a chance
at doing real science. The rest is just MSNBC graphics. :-)

Cheers,

David

--
David Fanning, Ph.D.
Fanning Software Consulting
Phone: 970-221-0438 E-Mail: davidf@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Toll-Free IDL Book Orders: 1-888-461-0155
  Switch to threaded view of this topic Create a new topic Submit Reply
Previous Topic: Re: IDL compatibility
Next Topic: Tickmarks and ticklabels

-=] Back to Top [=-
[ Syndicate this forum (XML) ] [ RSS ] [ PDF ]

Current Time: Thu Oct 09 06:43:08 PDT 2025

Total time taken to generate the page: 0.48012 seconds