Folks,

I have long thought that the IDL gridding routine, GridData,

to be one of IDL's most powerful and useful routines. Perhaps

taking its place among the likes of Histogram and Value_Locate.

Well, it *would* be powerful and useful if I could ever

get the damn thing to work. But, alas, I never have been

able to accomplish this simple feat.

I've decided to come clean about my abysmal failure

and ask for your help.

I ran into the perfect test case this week. A simple nearest

neighbor gridding problem that I know how to solve in two

completely independent ways, each producing identical

results. I *know* what I am doing here and I am

*supremely* confident in the results. "And," I thought,

"it is so simple, I could do this in GridData!"

Not. :-(

I've explained the problem and put some data here on

my web page:

http://www.yqcomputer.com/

I would be *extremely* grateful to anyone who can take

me by the hand and lead me to the promised land.

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

n Apr 17, 7:19m, David Fanning < XXXX@XXXXX.COM > wrote:

Hi David,

I've been trying something very similar recently. I think the

confusion is with the map_project operation - what you want is to

convert the UV coordinates of the stereo projection into lat,lon

values, and not the lat,lon into equirectangular UV. In the end,

though, you don't need Griddata gor it - hope that's not a

disappointment!

You can recreate the map you want (http://hrscview.fu-berlin.de/mex4/

software/other/out.png - sorry about the shocking colour!) like this:

pro tmp_fanning_map

im=fltarr(144,73)

openr,1,"D:\mydocs\work\2010-04-18 fanning\usegriddata.dat"

readu,1,im

close,1

im=reverse(im,2)

lat=gm_scl(indgen(73),out_range=[-90.,90.])

lon=gm_scl(indgen(144),out_range=[0.,360])

map=map_proj_init('Stereographic', center_lon=-45, center_lat=90,

sphere_radius=6378273.00)

sz=[304,448]

xr=[-385,375]*1e4

yr=[-535,585]*1e4

x=gm_scl(indgen(sz[0]),out_range=xr)

y=gm_scl(indgen(sz[1]),out_range=yr)

q=lindgen(product(sz))

qx=q mod sz[0]

qy=q/sz[0]

lonlat=map_proj_inverse(x[qx],y[qy],map_structure=map)

lon0=reform(lonlat[0,*])

lat0=reform(lonlat[1,*])

x0=wrap360(lon0)/360.*144

y0=(lat0[q]+90.)/180.*73.

out=fltarr(sz)

out[q]=im[x0[q],y0[q]]

device,decomposed=0

loadct,11

tvscl,out

end

You'll need these, too:

function gm_scl,x,in_range=in_range,out_range=out_range

;more powerful bytscl function - type taken from out_range if

present, otherwise x

;in_range - range of input data to be stretched

;out_range - output range

tname=size(keyword_set(out_range)?out_range:x,/tname)

type=size(keyword_set(out_range)?out_range:x,/type)

mn=min(x,max=mx)

if keyword_set(in_range) then begin

mn=in_range[0]

mx=in_range[1]

endif

if ~keyword_set(out_range) then begin

eps=1d-9

case type of

"UINT":out_range=[0,65536-eps]

"INT":out_range=[-32768,32768-eps]

"BYTE":out_range=[0,256-eps]

else:out_range=[0,100-eps]

endcase

endif

y=(double(x)-mn)/(mx-mn)

y=>>0<<1d

out=out_range[0]+y*(out_range[1]-out_range[0])

return,fix(out,type=type)

end

function wrap360,a ;make -179 and +179 close neighbours

return,(a+360) mod 360

end

Hi David,

I've been trying something very similar recently. I think the

confusion is with the map_project operation - what you want is to

convert the UV coordinates of the stereo projection into lat,lon

values, and not the lat,lon into equirectangular UV. In the end,

though, you don't need Griddata gor it - hope that's not a

disappointment!

You can recreate the map you want (http://hrscview.fu-berlin.de/mex4/

software/other/out.png - sorry about the shocking colour!) like this:

pro tmp_fanning_map

im=fltarr(144,73)

openr,1,"D:\mydocs\work\2010-04-18 fanning\usegriddata.dat"

readu,1,im

close,1

im=reverse(im,2)

lat=gm_scl(indgen(73),out_range=[-90.,90.])

lon=gm_scl(indgen(144),out_range=[0.,360])

map=map_proj_init('Stereographic', center_lon=-45, center_lat=90,

sphere_radius=6378273.00)

sz=[304,448]

xr=[-385,375]*1e4

yr=[-535,585]*1e4

x=gm_scl(indgen(sz[0]),out_range=xr)

y=gm_scl(indgen(sz[1]),out_range=yr)

q=lindgen(product(sz))

qx=q mod sz[0]

qy=q/sz[0]

lonlat=map_proj_inverse(x[qx],y[qy],map_structure=map)

lon0=reform(lonlat[0,*])

lat0=reform(lonlat[1,*])

x0=wrap360(lon0)/360.*144

y0=(lat0[q]+90.)/180.*73.

out=fltarr(sz)

out[q]=im[x0[q],y0[q]]

device,decomposed=0

loadct,11

tvscl,out

end

You'll need these, too:

function gm_scl,x,in_range=in_range,out_range=out_range

;more powerful bytscl function - type taken from out_range if

present, otherwise x

;in_range - range of input data to be stretched

;out_range - output range

tname=size(keyword_set(out_range)?out_range:x,/tname)

type=size(keyword_set(out_range)?out_range:x,/type)

mn=min(x,max=mx)

if keyword_set(in_range) then begin

mn=in_range[0]

mx=in_range[1]

endif

if ~keyword_set(out_range) then begin

eps=1d-9

case type of

"UINT":out_range=[0,65536-eps]

"INT":out_range=[-32768,32768-eps]

"BYTE":out_range=[0,256-eps]

else:out_range=[0,100-eps]

endcase

endif

y=(double(x)-mn)/(mx-mn)

y=>>0<<1d

out=out_range[0]+y*(out_range[1]-out_range[0])

return,fix(out,type=type)

end

function wrap360,a ;make -179 and +179 close neighbours

return,(a+360) mod 360

end

Greg writes:

OK, thanks. Now I know *three* ways to get the right

answer, but I still don't know how to use GridData.

Surely, *someone* has figured out how to use this

program!

Cheers,

David

P.S. I know I'm too circumspect at times, so let me

just say again, this is not a question about getting

the right answer, I know how to do that. This is a

question about getting GridData to give me the right

answer, which I presume is its purpose in life. :-)

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

OK, thanks. Now I know *three* ways to get the right

answer, but I still don't know how to use GridData.

Surely, *someone* has figured out how to use this

program!

Cheers,

David

P.S. I know I'm too circumspect at times, so let me

just say again, this is not a question about getting

the right answer, I know how to do that. This is a

question about getting GridData to give me the right

answer, which I presume is its purpose in life. :-)

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

Hi David, I have no problems with GRIDDATA; take a look at the code.

The only problem I had was the triangulate function - you might have

problems with collinear points on the poles if you don't remove

them).

Cheers, Klemen

pro tmp_fanning_map

;input size in x and y direction

sx = 144

sy = 73

;read input data

im=fltarr(sx,sy)

openr,1,"usegriddata.dat"

readu,1,im

close,1

im=reverse(im,2)

;generate input lon and lat array

im_lat = rebin(reform(findgen(sy)/(sy-1)*180.-90., 1, sy), sx, sy)

im_lon = rebin(findgen(sx)/sx*360., sx, sy)

;reduce the dimension in y directon (otherwise problems with colinear

points on the poles)

sy = sy - 2

im_lon = im_lon[*,1:sy]

im_lat = im_lat[*,1:sy]

im = im[*,1:sy]

;Polar projection on WGS84

map=map_proj_init(106, DATUM=8, /GCTP, center_lon=-45.,

center_lat=90.)

;transform input coorduiante arrays into vector

v_x = transpose(im_lon[*])

v_y = transpose(im_lat[*])

point_prj = MAP_PROJ_FORWARD(v_x, v_y, MAP_STRUCTURE=map)

;Make triangles

TRIANGULATE, point_prj[0,*], point_prj[1,*], Trng, TOLERANCE=1.

im_prj = GRIDDATA(point_prj[0,*], point_prj[1,*], im[*], $

/NEAREST_NEIGHBOR, DELTA=[25000.,25000.], TRIANGLES=Trng, $

DIMENSION=[304,448], START=[-3850000., -5350000.])

im_prj = reverse(im_prj, 2)

save, im_prj, file='faning.sav'

device,decomposed=0

loadct,11

tvscl,im_prj

end

The only problem I had was the triangulate function - you might have

problems with collinear points on the poles if you don't remove

them).

Cheers, Klemen

pro tmp_fanning_map

;input size in x and y direction

sx = 144

sy = 73

;read input data

im=fltarr(sx,sy)

openr,1,"usegriddata.dat"

readu,1,im

close,1

im=reverse(im,2)

;generate input lon and lat array

im_lat = rebin(reform(findgen(sy)/(sy-1)*180.-90., 1, sy), sx, sy)

im_lon = rebin(findgen(sx)/sx*360., sx, sy)

;reduce the dimension in y directon (otherwise problems with colinear

points on the poles)

sy = sy - 2

im_lon = im_lon[*,1:sy]

im_lat = im_lat[*,1:sy]

im = im[*,1:sy]

;Polar projection on WGS84

map=map_proj_init(106, DATUM=8, /GCTP, center_lon=-45.,

center_lat=90.)

;transform input coorduiante arrays into vector

v_x = transpose(im_lon[*])

v_y = transpose(im_lat[*])

point_prj = MAP_PROJ_FORWARD(v_x, v_y, MAP_STRUCTURE=map)

;Make triangles

TRIANGULATE, point_prj[0,*], point_prj[1,*], Trng, TOLERANCE=1.

im_prj = GRIDDATA(point_prj[0,*], point_prj[1,*], im[*], $

/NEAREST_NEIGHBOR, DELTA=[25000.,25000.], TRIANGLES=Trng, $

DIMENSION=[304,448], START=[-3850000., -5350000.])

im_prj = reverse(im_prj, 2)

save, im_prj, file='faning.sav'

device,decomposed=0

loadct,11

tvscl,im_prj

end

Klemen writes:

Ah, ha! Thanks, Klemen.

This is similar to what I had originally done, but I couldn't

overcome that darn collinear points problem. Eventually, I made

triangles from the original data XY and that "worked", so I thought

that must be right. It never occurred to me that the pole points

could be the problem.

Just one correction to your code. You set up a polar

stereographic map projection like this:

map=map_proj_init(106, DATUM=8, /GCTP, center_lon=-45., center_lat=90.)

This is incorrect. In the polar stereographic projection (only!!)

the CENTER_LAT keyword is actually the TRUE_SCALE_LATITUDE

keyword. (This is being fixed, thank goodness, in IDL 8.0, since

*everyone* makes this mistake.)

http://www.yqcomputer.com/

The projection should really be set up like this:

map=map_proj_init(106, DATUM=8, /GCTP, center_lon=-45., center_lat=70.)

I'm really grateful for your help. :-)

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

Ah, ha! Thanks, Klemen.

This is similar to what I had originally done, but I couldn't

overcome that darn collinear points problem. Eventually, I made

triangles from the original data XY and that "worked", so I thought

that must be right. It never occurred to me that the pole points

could be the problem.

Just one correction to your code. You set up a polar

stereographic map projection like this:

map=map_proj_init(106, DATUM=8, /GCTP, center_lon=-45., center_lat=90.)

This is incorrect. In the polar stereographic projection (only!!)

the CENTER_LAT keyword is actually the TRUE_SCALE_LATITUDE

keyword. (This is being fixed, thank goodness, in IDL 8.0, since

*everyone* makes this mistake.)

http://www.yqcomputer.com/

The projection should really be set up like this:

map=map_proj_init(106, DATUM=8, /GCTP, center_lon=-45., center_lat=70.)

I'm really grateful for your help. :-)

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

Hi David,

thank you for TRUE_SCALE_LATITUDE tip, I wasn't aware of it!

Cheers, Klemen

thank you for TRUE_SCALE_LATITUDE tip, I wasn't aware of it!

Cheers, Klemen

Klemen writes:

I've discovered a couple more interesting facts about

this process this morning. It turns out that it is

mostly Triangulate that is giving me problems. I've

found I do NOT have to exclude any values to produce

the proper triangulation, and that the "co-linear"

problem occurs on my Windows box, but not my Linux

box. On Windows, setting the TOLERANCE keyword to 1

appears to solve the problem.

Also interesting is that there is a small gap (mostly

camouflaged in my web article) where the longitude vector

wraps around from 257.5 to 0 degrees. This is especially

apparent in the filled contour method, and less apparent

in the NSIDC regrid method. It appears to disappear

completely in the GridData method, perhaps justifying

my confidence in its power, if we can just learn to

harness it. :-)

I'll update my web page article sometime soon. But

managing to do this with GridData opens up a path

I have been searching for for at least the last

two years!

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thue. ("Perhaps thos speakest truth.")

I've discovered a couple more interesting facts about

this process this morning. It turns out that it is

mostly Triangulate that is giving me problems. I've

found I do NOT have to exclude any values to produce

the proper triangulation, and that the "co-linear"

problem occurs on my Windows box, but not my Linux

box. On Windows, setting the TOLERANCE keyword to 1

appears to solve the problem.

Also interesting is that there is a small gap (mostly

camouflaged in my web article) where the longitude vector

wraps around from 257.5 to 0 degrees. This is especially

apparent in the filled contour method, and less apparent

in the NSIDC regrid method. It appears to disappear

completely in the GridData method, perhaps justifying

my confidence in its power, if we can just learn to

harness it. :-)

I'll update my web page article sometime soon. But

managing to do this with GridData opens up a path

I have been searching for for at least the last

two years!

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thue. ("Perhaps thos speakest truth.")

David Fanning writes:

Just one more surprise to report. If I use the GridData

Natural Neighbor method (instead of the Nearest Neighbor

method), the result is indistiguishable from the Contour

method!

I've updated my web page to report the new results:

http://www.yqcomputer.com/

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thue. ("Perhaps thos speakest truth.")

Just one more surprise to report. If I use the GridData

Natural Neighbor method (instead of the Nearest Neighbor

method), the result is indistiguishable from the Contour

method!

I've updated my web page to report the new results:

http://www.yqcomputer.com/

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thue. ("Perhaps thos speakest truth.")

David Fanning writes:

Having GridData working is like living in the land of

milk and honey. But there is one disturbing fact that

lurks like a snake in the garden. I get different

results from GridData depending upon what machine

I run the identical code on. :-(

In particular, the Modifed Shepard's algorithm blows

up (no other way to describe it) on my Windows 64-bit

machine, while it seems to work reasonably well on

my 32-bit LINUX machine.

Can anyone think of any reasonable explanation for this?

I show some results at the end of this article:

http://www.yqcomputer.com/

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

Having GridData working is like living in the land of

milk and honey. But there is one disturbing fact that

lurks like a snake in the garden. I get different

results from GridData depending upon what machine

I run the identical code on. :-(

In particular, the Modifed Shepard's algorithm blows

up (no other way to describe it) on my Windows 64-bit

machine, while it seems to work reasonably well on

my 32-bit LINUX machine.

Can anyone think of any reasonable explanation for this?

I show some results at the end of this article:

http://www.yqcomputer.com/

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

> > ttp://www.dfanning.com/code_tips/usegriddata.html >> >> Having GridData working is like living in the land of >> milk and honey. But there is one disturbing fact that >> lurks like a snake in the garden. I get different >> results from GridData depending upon what machine >> I run the identical code on. :-( >> >> In particular, the Modifed Shepard's algorithm blows >> up (no other way to describe it) on my Windows 64-bit >> machine, while it seems to work reasonably well on >> my 32-bit LINUX machine. >> >> Can anyone think of any reasonable explanation for this? >> I show some results at the end of this article: >> >> ttp://www.dfanning.com/code_tips/usegriddata.html> >

Hi David,

It doesn't look to me like you are running the coordinates through

GRID_INPUT before passing them to GRIDDATA. In the past I found that

(generally) removed all subsequent potholes from the GRIDDATA process.

Cheers,

Ben

ben.bighair writes:

Well, I tried this:

Grid_Input, x, y, air, x1, y1, air1

Triangulate, x1, y1, triangles1, TOLERANCE=1.0

griddedData = GRIDDATA(x1, y1, air1, $

/SHEPARDS, DELTA=[25000.,25000.], TRIANGLES=triangles1, $

DIMENSION=[304,448], START=[-3850000., -5350000.])

Now the Shepard's plot "blows up", but in a more orderly

fashion, with straight lines. :-)

Still works on my LINUX box, though, which I find most

disconcerting.

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

Well, I tried this:

Grid_Input, x, y, air, x1, y1, air1

Triangulate, x1, y1, triangles1, TOLERANCE=1.0

griddedData = GRIDDATA(x1, y1, air1, $

/SHEPARDS, DELTA=[25000.,25000.], TRIANGLES=triangles1, $

DIMENSION=[304,448], START=[-3850000., -5350000.])

Now the Shepard's plot "blows up", but in a more orderly

fashion, with straight lines. :-)

Still works on my LINUX box, though, which I find most

disconcerting.

Cheers,

David

--

David Fanning, Ph.D.

Fanning Software Consulting, Inc.

Coyote's Guide to IDL Programming: http://www.yqcomputer.com/

Sepore ma de ni thui. ("Perhaps thou speakest truth.")

2. griddata 'v4' interpolation

3. griddata produces NaN ? what's wrong

4. griddata with concave regions

5. what's the diffirent from griddata and interp2

6. griddata error: qhull precision error, ERRONEOUS FACET

9. how to get values of arbitrary positions after griddata

10. griddata help: want max in cell instead of linear interp

11. Using griddata with method 'v4'

12. SAR polar formatting: using griddata

13. griddata

14. Using GRIDDATA on a closed surface

15. griddata function () in matlab

11 post • Page:**1** of **1**