| Re: Overplot nice looking globe on 2d satellite images [message #52293] |
Mon, 29 January 2007 09:31 |
Brian Larsen
Messages: 270
Registered: June 2006
|
Senior Member |
|
|
Oops, I found an error that is pretty important.
The following two lines make the Earth turn the wrong way
> rot = (rot[2]+rot[3]/60.+rot[4]/3600.)/24.*360 - 270
> ENDIF ELSE rot = time/24. * 360 - 270
replace them with:
rot = -(rot[2]+rot[3]/60.+rot[4]/3600.)/24.*360 + 90
ENDIF ELSE rot = -time/24. * 360 + 90
for an Earth that turns the correct direction.
Brian
------------------------------------------------------------ ---------
Brian A. Larsen
Dept. of Physics
Space Science and Engineering Lab (SSEL)
Montana State University - Bozeman
Bozeman, MT 59717
|
|
|
|
| Re: Overplot nice looking globe on 2d satellite images [message #52294 is a reply to message #52293] |
Mon, 29 January 2007 09:02  |
James Kuyper
Messages: 425
Registered: March 2000
|
Senior Member |
|
|
David Fanning wrote:
> Brian Larsen writes:
>
>> I imagine this is not too hard but I haven't been able to figure it
>> out. What I want to do is in the center of a polar contour plot
>> overplot a nice looking globe that can be turned and tilted to the
>> correct perspective and shaded to show where the sun is.
>
> In my experience, any map projection that overlaps the pole
> and includes the sun is trouble:
>
> http://www.dfanning.com/map_tips/terminator.html
>
> But maybe this is because I don't understand spherical
> geometry worth a damn. :-)
I've just taken a quick look at that code. My initial impression is
that you shouldn't have a problem if you check whether a give point is
inside or outside the terminator polygon, by performing the check in
projection coordinates, rather than latitude/longitude. However, that
depends upon the details of the Inside() subroutine, about which I
know nothing that isn't implied by it's name.
Another, more direct approach is possible. Convert both the ground
location and the sub-solar location into unit vectors in earth-
centered rotating coordinates (if you want to take into consideration
the fact that the earth isn't perfectly spherical, the equations get
more complicated):
<cos(lat)*cos(lon), cos(lat)*sin(lon), sin(lat)>
There are many different definitions of the day/night distinction,
depending upon whether you want nightfall to occur when the center of
the sun is on the horizon, or when the entire sun is below the
horizon. It also depends upon whether you want to use the geometric
position of the sun, or the shifted position due to atmospheric
refraction. All of these issues can be covered by choosing a suitable
value for the terminator_angle, the angle between the sun's position
and the terminator; if you're using the center of the sun with no
correction for refraction, the terminator_angle should be pi/2
radians, or 90 degrees..
If the dot product of the ground position unit vector and the sub-
solar point's unit vector is greater than cos(terminator_angle), then
the ground position is on the day side of the terminator. If it's
less, then it's on the night side.
|
|
|
|
| Re: Overplot nice looking globe on 2d satellite images [message #52295 is a reply to message #52294] |
Mon, 29 January 2007 07:58 |
Brian Larsen
Messages: 270
Registered: June 2006
|
Senior Member |
|
|
Como la mayor parte de nosotros mi español es muy malo y se retarda,
pero te pareces gozar de tu viaje a España
:)
Brian
------------------------------------------------------------ ---------
Brian A. Larsen
Dept. of Physics
Space Science and Engineering Lab (SSEL)
Montana State University - Bozeman
Bozeman, MT 59717
On Jan 28, 5:39 pm, David Fanning <d...@dfanning.com> wrote:
> Brian Larsen writes:
>> At the risk of it being hacked apart, to finish off the thread for any
>> to come later I will attach the final code that I ended up using as a
>> procedure plot_earth.proAh, bueno! Verdad?
>
> Saludo,
>
> David
> --
> David Fanning, Ph.D.
> Fanning Software Consulting, Inc.
> Coyote's Guide to IDL Programming:http://www.dfanning.com/
|
|
|
|
|
|
| Re: Overplot nice looking globe on 2d satellite images [message #52304 is a reply to message #52303] |
Sun, 28 January 2007 14:46 |
Brian Larsen
Messages: 270
Registered: June 2006
|
Senior Member |
|
|
In my best Homer Simpson, "Doh!!"
At the risk of it being hacked apart, to finish off the thread for any
to come later I will attach the final code that I ended up using as a
procedure plot_earth.pro
;+
; NAME:
; plot_earth
;
;
; PURPOSE:
; overplot the earth onto a plot measured in Re
;
;
; INPUTS:
; none
;
;
; OPTIONAL INPUTS:
; none
;
;
; KEYWORD PARAMETERS:
; FANCY - plot the earth as a globe
; BLANK - for non fancy output blank out the region, ingnored when
; fancy is specified
; TIME - the time of the fancy earth to plot, defults to zero
; This can be specifed in fractional hours (e.g. 12.345) or as
; an EUVtime (e.g. 20012001234 yyyydddhhmmss)
;
;
; OUTPUTS:
; the earth onto the current plot
;
;
; OPTIONAL OUTPUTS:
; none
;
;
; COMMON BLOCKS:
; none
;
;
; SIDE EFFECTS:
; - puts the earth on the plot
; - messed up the colortable, need to reload ct after a run of this
;
; RESTRICTIONS:
; - Only works well for iso plots measured in Re from the earth.
Could
; be modified to work in meters or the like at a later date
; - Only works for directly over the North Pole as that I what I am
doing
; now
;
;
; EXAMPLE:
; loadct, 13
; contour, dist(10), findgen(10)-5, findgen(10)-5, /fill, nlevels=30, /
iso
; plot_earth, /fancy, time=12.345
;
;
; MODIFICATION HISTORY:
;
; Sun Jan 28 15:41:46 2007, Brian Larsen
; <larsen@ssel.montana.edu>
;
; written and tested with help from comp.lang.idl-pvwave
; http://groups.google.com/group/comp.lang.idl-pvwave/
browse_thread/thread/4ba6a41393fc8f6f/?hl=en#
;
;-
pro plot_earth, BLANK=blank, FANCY=fancy, TIME=time
IF NOT KEYWORD_SET(fancy) THEN BEGIN
IF KEYWORD_SET(blank) THEN BEGIN
circ_r = fltarr(200)
circ_r[*] = 1
circ = findgen(200)*2*!pi/200
x_=circ_r * sin(circ)
y_=circ_r * cos(circ)
polyfill, x_, y_, color=0
ENDIF
rad_ = 2.*!pi*findgen(100)/100.
earth_ = fltarr(100)
earth_[*] = 1
oplot, /polar, earth_[*], 2.*rad_[*]
plots, [0,1],[0,0], linestyle=2
ENDIF ELSE BEGIN
;; for time we are expecting either an euv time or factional hours
;; - an euvtime will be a string
IF N_ELEMENTS(time) EQ 0 THEN rot = 0 ELSE BEGIN
IF size(time, /type) EQ 7 THEN BEGIN
rot = euv_date2arr( time)
rot = (rot[2]+rot[3]/60.+rot[4]/3600.)/24.*360 - 270
ENDIF ELSE rot = time/24. * 360 - 270
ENDELSE
;; day side
pos1 = convert_coord([0,1], [-1,1], /data, /to_normal)
;; night side
pos2 = convert_coord([-1,0], [-1,1], /data, /to_normal)
;; day side
loadct, 12
MAP_SET,90,0,rot,/ORTHOGRAPHIC,/ISOTROPIC, /CONTINENTS,/HORIZON, $
E_continents={FILL:1, color:23}, $
position=[pos1[0,0], pos1[1,0], pos1[0,1], pos1[1,1]], $
/noerase, /noborder, $
E_HORIZON={FILL:1, COLOR:100}, $
limit=[0,rot,90,rot+180]
;; night side
loadct, 0
MAP_SET,90,0,rot,/ORTHOGRAPHIC,/ISOTROPIC, /CONTINENTS,/HORIZON, $
E_continents={FILL:1, color:119}, $
position=[pos2[0,0], pos2[1,0], pos2[0,1], pos2[1,1]], $
/noerase, /noborder, $
E_HORIZON={FILL:1, COLOR:33}, $
limit=[0,rot+180,90,rot]
ENDELSE
END
------------------------------------------------------------ ---------
Brian A. Larsen
Dept. of Physics
Space Science and Engineering Lab (SSEL)
Montana State University - Bozeman
Bozeman, MT 59717
|
|
|
|
| Re: Overplot nice looking globe on 2d satellite images [message #52305 is a reply to message #52304] |
Sun, 28 January 2007 13:42 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Brian Larsen writes:
> Well adding more to my own posts, interestingly the map is not in the
> right spot. Any insight here?
>
> Looking at the plot generated (sorry for the terrible colors) it is
> pretty clear that the earth is not centered at 0,0 is there some trick
> I don't know here?
>
> loadct, 11
> contour, dist(10), findgen(10)-5, findgen(10)-5, /fill, nlevels=30, /
> iso
> pos = convert_coord([-1,1], [-1,1], /data, /to_normal)
> loadct, 12
> MAP_SET,90,0,0,/ORTHOGRAPHIC,/ISOTROPIC, /CONTINENTS,/HORIZON, $
> E_continents={FILL:1, color:23}, $
> position=[pos[0,0], pos[0,1], pos[1,0], pos[1,1]], $
> /noerase, /noborder, $
> E_HORIZON={FILL:1, COLOR:100}
You set up your position wrong. It should be set like this:
position=[pos[0,0], pos[1,0], pos[0,1], pos[1,1]], $
Then, it is right where it is suppose to be. :-)
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
|
|
|
| Re: Overplot nice looking globe on 2d satellite images [message #52306 is a reply to message #52305] |
Sun, 28 January 2007 10:55 |
Brian Larsen
Messages: 270
Registered: June 2006
|
Senior Member |
|
|
Well adding more to my own posts, interestingly the map is not in the
right spot. Any insight here?
Looking at the plot generated (sorry for the terrible colors) it is
pretty clear that the earth is not centered at 0,0 is there some trick
I don't know here?
loadct, 11
contour, dist(10), findgen(10)-5, findgen(10)-5, /fill, nlevels=30, /
iso
pos = convert_coord([-1,1], [-1,1], /data, /to_normal)
loadct, 12
MAP_SET,90,0,0,/ORTHOGRAPHIC,/ISOTROPIC, /CONTINENTS,/HORIZON, $
E_continents={FILL:1, color:23}, $
position=[pos[0,0], pos[0,1], pos[1,0], pos[1,1]], $
/noerase, /noborder, $
E_HORIZON={FILL:1, COLOR:100}
Brian
------------------------------------------------------------ ---------
Brian A. Larsen
Dept. of Physics
Space Science and Engineering Lab (SSEL)
Montana State University - Bozeman
Bozeman, MT 59717
|
|
|
|
| Re: Overplot nice looking globe on 2d satellite images [message #52307 is a reply to message #52306] |
Sun, 28 January 2007 10:32 |
Brian Larsen
Messages: 270
Registered: June 2006
|
Senior Member |
|
|
Thanks, this is huge forward progress
I notice that the output from convert_coord doesn't match position
keyword (typical idl annoyance)
something like this is what you had in mind? (looking straight down
from over the pole and w/o worrying about shading yet)
IDL> plot, findgen(11)-5, findgen(11)-5, /nodata, /iso
IDL> print, convert_coord([-1,1], [-1,1], /data, /to_normal)
0.356255 0.441672 0.00000
0.468755 0.591672 0.00000
IDL> MAP_SET,90,0,0,/ORTHOGRAPHIC,/ISOTROPIC, /CONTINENTS,/HORIZON, $
E_continents={FILL:1},
position=[0.356255,0.441672,0.468755,0.591672], /noerase, /noborder
(I will of course code in the position as variables since it will
change)
thanks much
Brian
------------------------------------------------------------ ---------
Brian A. Larsen
Dept. of Physics
Space Science and Engineering Lab (SSEL)
Montana State University - Bozeman
Bozeman, MT 59717
|
|
|
|
| Re: Overplot nice looking globe on 2d satellite images [message #52308 is a reply to message #52307] |
Sun, 28 January 2007 07:47 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Brian Larsen writes:
> I imagine this is not too hard but I haven't been able to figure it
> out. What I want to do is in the center of a polar contour plot
> overplot a nice looking globe that can be turned and tilted to the
> correct perspective and shaded to show where the sun is.
In my experience, any map projection that overlaps the pole
and includes the sun is trouble:
http://www.dfanning.com/map_tips/terminator.html
But maybe this is because I don't understand spherical
geometry worth a damn. :-)
> One thing that seems hard to do with the map routines (that I am a
> major novice at using as I really only look at space data with a
> little earth in the center) are that I see no easy way to make a globe
> take up -1,1 in a plot and leave the other stuff plotted around it.
Humm. Have you tried the keyword POSITION, NOERASE, and NOMARGIN?
It seems to me that would be about as easy as it comes. You might
need to use CONVERT_COORD to convert data coordinate space into
normalized space for the POSITION keyword, but that is the only
complication I can think of.
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
