Grady Daub <gadZOOKS8371@garnet.acns.fsummer.edu> wrote:
> Martin Schultz wrote:
> I've tried TRIANGULATE and TRIGRID and SPH_SCAT (the latter exactly as shown on dfanning.com, except with my own data)
> and the results are weird. Is SPH_SCAT, when applied to data with max/min of around 200/-200, supposed to produce
> output past 50000? That's the weird part.
I have had problems with sph_scat, especially near the (South) pole - see my recent
post with sph_scat in the title and the reply I received which might help you.
For my part, to get acceptable interpolation near the poles, I had to go via device
coordinates. The following code fragment might help (or not...)
;
; Use sph_scat
;
if (keyword_set(sphere)) then begin
lo=data(0,*)
la=data(1,*)
d2=data(2,*)
res=sph_scat(lo,la,d2,gs=gs,bounds=[xr(0),yr(0),xr(1),yr(1)] )
;
; Go via device coords
;
endif else if (keyword_set(by_dev)) then begin
; Convert from data coords to device coords. Its up to you to make sure that the
; mapping you want has been set up - eg by set_ps_map
xy=convert_coord(data(0,*),data(1,*),/data,/to_dev)
triangulate,xy(0,*),xy(1,*),triangles
; Now, regrid onto a grid regular in device space
res1=trigrid(xy(0,*),xy(1,*),data(2,*),triangles,nx=int_res, ny=int_res,xg=xg,yg=yg)
; Now, generate the lat-lon grid we want and convert it too into device coord
x=xc(pp,/arr,/as) & y=yc(pp,/arr)
xy1=convert_coord(x,y,/data,/to_dev)
; Now, interpolate from our new regular grid onto the xformed lat-lon grid...
; Since interpolate assumes that the input grid has coords 0..n-1 we need to xform
; the coords into this range...
; Compute scaling that *would* xform xg, yg to 0...int_res-1
xr=makerange(xg) & yr=makerange(yg)
mn=[xr(0),yr(0)] & sc=1./[xr(1)-xr(0),yr(1)-yr(0)]
xy1a=xy1 & for i=0,1 do xy1a(i,*)=(int_res-1)*(xy1(i,*)-mn(i))*sc(i)
res=interpolate(res1,xy1a(0,*),xy1a(1,*))
res=reform(res,pp.lbnpt,pp.lbrow)
;
; Or just do it in lat-lon coords - OK if not too near the pole
;
endif else begin
triangulate,data(0,*),data(1,*),triangles
res=trigrid(data(0,*),data(1,*),data(2,*),triangles,gs,[xr(0 ),yr(0),xr(1),yr(1)])
endelse
>> Finally (to complete the most prominent issues with contours on maps),
>> if you plan to plot filled contours, there are occasions when you need
>> the /CELL_FILL keyword ...
> Where is /CELL_FILL documented? It's not it the index and I've not found it anywhere near the CONTOUR section.
/cell_fill was supposed to be obsolete and was actually removed in one release
(5.0?) - it fills cell-by-cell not contour-by-contour. I too find that it needs
to be used sometimes, which is annoying because it is slow and doesn't work with
patterns. RSI please note and fix this!
- William
--
William M Connolley | wmc@bas.ac.uk | http://www.nbs.ac.uk/public/icd/wmc/
Climate Modeller, British Antarctic Survey | Disclaimer: I speak for myself
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
