| Map Projections and Contour Plots. [message #15827] |
Thu, 10 June 1999 00:00  |
Grady Daub
Messages: 22
Registered: June 1999
|
Junior Member |
|
|
I have data that corresponds to latitude and longitude.
example. data=40 , lat=20, lon=20
How would I...
...make the correct arrays expected by CONTOUR?
/IRREGULAR and TRIANGULATION work, but, I think they're causing
problems.
...smooth the contours?
The documentation says something about MIN_CURVE_SURF, but, it seems to
need the same format array as CONTOUR.
...plot the contours correctly over a map projection?
When I try to CONTOUR over a map projection, I get some lines across the
plot. The lines look like the contour tried to plot from, say, longitude
180 to -180.
So far, I'm only able to put a map ONTO a CONTOUR, not the other way
around. This doesn't help me when I want to use, say, /MOLLWEIDE.
-Grady Daub
GFDI
To reply be email, please, remove "ZOOKS" and "MMER".
|
|
|
|
| Re: Map Projections and Contour Plots. [message #15832 is a reply to message #15827] |
Wed, 16 June 1999 00:00  |
Martin Schultz
Messages: 515
Registered: August 1997
|
Senior Member |
|
|
Grady Daub wrote:
>
> Martin Schultz wrote:
>
> Your convert_lon program was not attached.
sorry! this time it should be ...
>
> Would your program make my non-monotonic lat,lon,data arrays in a form expected by CONTOUR ?
no. it just converts lon value by lon value to fall into the range of
-180 to +180 or 0 to 360. You woul dneed to use the SORT function to get
monotonic values first.
>
> 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 haven't used triangulate and the likes for a long time, and I must
have deleted the code where I last did it. But I remember it had helped
me a great deal to actually plot the triangles returned by triangulate
in order to believe what IDL (and I) were doing...
(The CELL_FILL question was already answered by wmc)
Regards,
Martin
--
|||||||||||||||\\\\\\\\\\\\\-------------------///////////// //|||||||||||||||
Martin Schultz, DEAS, Harvard University, 29 Oxford St., Pierce 109,
Cambridge, MA 02138 phone (617) 496 8318 fax (617) 495 4551
e-mail mgs@io.harvard.edu web http://www-as/people/staff/mgs/
; $Id: convert_lon.pro,v 1.11 1999/05/20 00:58:48 mgs Exp $
;----------------------------------------------------------- --
;+
; NAME:
; CONVERT_LON
;
; PURPOSE:
; Convert longitudes from -180..180 to 0..360 or vice
; versa.
;
; CATEGORY:
; Tools
;
; CALLING SEQUENCE:
; CONVERT_LON,data,names,Pacific=Pacific,Atlantic=Atlantic, $
; minval=minval
;
; INPUTS:
; DATA -> A data array (lines,vars) or vector containing
; longitude data. If DATA is a 2D array, the NAMES
; parameter must be given to identify the LONgitude variable.
;
; NAMES -> A string list of variable names. The longitude data
; must be labeled 'LON', unless specified with the LONNAME
; keyword. The NAMES parameter is not needed, if a data
; vector is passed.
;
; KEYWORD PARAMETERS:
; PACIFIC -> Convert longitudes from -180..180 to 0..360
;
; ATLANTIC -> Convert from 0..360 to -180..180
;
; LONNAME -> Name of the longitude variable if a name other
; than 'LON' is used.
;
; OUTPUTS:
; The longitude column in the data array will be changed.
;
; SUBROUTINES:
;
; REQUIREMENTS:
;
; NOTES:
;
; EXAMPLE:
; londat = [ -180.,-179.,-0.1,0.1,179.,180.,270.,359.]
; CONVERT_LON,londat,/Pacific
; print,londat
;
; CONVERT_LON,londat,/Atlantic
; print,londat
;
; MODIFICATION HISTORY:
; mgs, 25 Aug 1998: VERSION 1.00
; mgs, 19 May 1999: - now makes sure that longitude range does
; not exceed -180..180 or 0..360
;
;-
; Copyright (C) 1998, Martin Schultz, Harvard University
; This software is provided as is without any warranty
; whatsoever. It may be freely used, copied or distributed
; for non-commercial purposes. This copyright notice must be
; kept with any copy of this software. If this software shall
; be used commercially or sold as part of a larger package,
; please contact the author to arrange payment.
; Bugs and comments should be directed to mgs@io.harvard.edu
; with subject "IDL routine convert_lon"
;----------------------------------------------------------- --
pro convert_lon,data,names,pacific=pacific,atlantic=atlantic, $
lonname=lonname
minval = -180.0001
if (n_elements(lonname) eq 0) then lonname = 'LON'
if (n_elements(data) lt 2) then return
; get size information of data and find LON column
s = size(data)
if (s[0] eq 1) then ind = 0 $ ; data is vector
else begin
; Find LON variable
ind = where(strupcase(names) eq lonname)
if (ind[0] lt 0) then begin
print,'*** CONVERT_LON: Cannot find ',lonname,' in data set!'
return
endif
endelse
; Atlantic: Convert longitudes greater 180 by subtracting 360
; also add N*360 to longitude values less than -180
if (keyword_set(Atlantic)) then begin
repeat begin
lon = data[*,ind[0]]
index = where(lon gt 180.,count)
if (index[0] ge 0) then data[index,ind[0]] = lon[index]-360.
endrep until(count eq 0)
repeat begin
lon = data[*,ind[0]]
index = where(lon lt -180.,count)
if (index[0] ge 0) then data[index,ind[0]] = lon[index]+360.
endrep until(count eq 0)
endif
; Pacific: convert negative longitudes by adding 360
; also subtract N*360 for longitude values greater than 360
if (keyword_set(Pacific)) then begin
repeat begin
lon = data[*,ind[0]]
; index = where(lon gt minval AND lon lt 0.)
index = where(lon lt 0., count)
if (index[0] ge 0) then data[index,ind[0]] = lon[index]+360.
endrep until(count eq 0)
repeat begin
lon = data[*,ind[0]]
index = where(lon gt 360., count)
if (index[0] ge 0) then data[index,ind[0]] = lon[index]-360.
endrep until(count eq 0)
endif
return
end
|
|
|
|
| Re: Map Projections and Contour Plots. [message #15841 is a reply to message #15827] |
Wed, 16 June 1999 00:00 |
wmc
Messages: 117
Registered: February 1995
|
Senior Member |
|
|
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
|
|
|
|
| Re: Map Projections and Contour Plots. [message #15858 is a reply to message #15827] |
Tue, 15 June 1999 00:00 |
Grady Daub
Messages: 22
Registered: June 1999
|
Junior Member |
|
|
Martin Schultz wrote:
> LONs and LATs have to be monotonically increasing for a contour plot.
> Default is -180 to 180 for a global map, but if you want to center your
> map say on the Pacific, you need to "convert" your longitude coordinates
> to e.g. 0 to 360. You can try my convert_lon program (attached) for
> this. It is called as
> convert_lon,lon [,/ATLANTIC] [,/PACIFIC]
> and makes sure you get longitudes modulo 360.
>
Your convert_lon program was not attached.
Would your program make my non-monotonic lat,lon,data arrays in a form expected by CONTOUR ?
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.
TRIANGULATE and TRIGRID just don't like me. :-(
What's the deal with FVALUE?
> 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 *AND* CONTOUR will use the first fill color for
> values ABOVE your first threshold, i.e. you should add a very low number
> to your C_LEVEL specification if you have one.
Where is /CELL_FILL documented? It's not it the index and I've not found it anywhere near the CONTOUR section.
-Grady Daub
(Remove MMER and ZOOKS to reply by email.)
|
|
|
|
| Re: Map projections [message #27968 is a reply to message #15827] |
Mon, 12 November 2001 08:16 |
Liam E. Gumley
Messages: 378
Registered: January 2000
|
Senior Member |
|
|
James Kuyper wrote:
> Is there any way to get at the coordinate conversion function for the
> currently active map projection? I'm talking about a function that takes
> a physical position in latitude and longitude coordinates, and converts
> it into an image position using either data, device, or normal
> coordinates. It doesn't matter which; once the position is in any one of
> those forms, I can get the others by using CONVERT_COORD. Such functions
> have to exist in order for IDL to perform mapping, but they don't seem
> to be publically exposed. I'd also like to have access to the inverse
> function, to get a lat/lon corresponding to an image position.
James,
When a map projection is active, longitude and latitude correspond to
x/y data coordinates. Therefore CONVERT_COORD will allow you to obtain
the device or normal coordinates for a given lon/lat. For example:
window, /free
map_set, 45, -90, scale=20e6, /usa
lon = -89.6
lat = 43.5
print, convert_coord(lon, lat, /data, /to_device)
326.465 217.622 0.000000
print, convert_coord(lon, lat, /data, /to_normal)
0.510102 0.425044 0.000000
Cheers,
Liam.
Practical IDL Programming
http://www.gumley.com/
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
