[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: warping continents to magnetic coordinates
Hi Ken
I've got some code that you may find useful. It
works in "geomagnetic dipole" coordinates as defined
in Russell or Hapgood. Fairly similar to CGM in
most places.
The basic idea is to pull data directly from the IDL
map files, then apply a coordinate transformation. The
rotation matrix is obtained from some *VERY* rough code
of mine called "geospace_convert". It seems to work
well for this purpose, but please use with caution.
Let me know if you have any trouble with this. It
*should* just compile and run. Again, this is not
high quality code, but hopefully you'll find it useful.
Brian
;1999/08/03 Brian Jackel
;
; This IDL script draws a continental map in dipole geomagnetic
coordinates.
; It takes a few key lines from "map_continents" and does coordinate
transformation
; using "geospace_convert".
;
name= FILEPATH('clow.ndx',SUBDIR=['resource','maps','low'])
OPENR,lun,fname,/XDR,/GET_LUN
nsegments=0L
READU,lun,nsegments
result= REPLICATE({fptr:0L, npts:0L, latmax:0.0, latmin:0.0,
lonmax:0.0, lonmin:0.0},nsegments)
READU,lun,result
FREE_LUN,lun
map_set,90,-50,79,/GNOMIC,LIMIT=[55,-120,85,190]
re= 6371.2
;
;
;rotation_matrix= FLTARR(3,3)
;rotation_matrix[0,0]= 1.0
;rotation_matrix[1,1]= 1.0
;rotation_matrix[2,2]= 1.0
CDF_EPOCH,Epoch,1990,1,1,23,0,0,/COMPUTE_EPOCH
dummy=
GEOSPACE_CONVERT([1,1,1],'geo','mag',ROTATION_MATRIX=rotation_matrix,
CDF_EPOCH_DATE=epoch)
fname= FILEPATH('clow.dat',SUBDIR=['resource','maps','low'])
OPENR,lun,fname,/XDR,/GET_LUN
FOR indx=0,nsegments-1 DO BEGIN
IF (result[indx].npts GE 2) THEN BEGIN
POINT_LUN,lun,result[indx].fptr
xy= FLTARR(2,result[indx].npts,/NOZERO)
READU,lun,xy
R= 1
lat_rad= REFORM(xy[0,*]) *!dtor
lon_rad= REFORM(xy[1,*]) *!dtor
temp= FLTARR(result[indx].npts,3)
temp(0,0)= R * cos(lat_rad) * cos(lon_rad)
temp(0,1)= R * cos(lat_rad) * sin(lon_rad)
temp(0,2)= R * sin(lat_rad)
temp= rotation_matrix##temp
l= SQRT( temp[*,0]^2 + temp[*,1]^2 )
xy[0,*]= ATAN(temp[*,2],l) * !radeg
xy[1,*]= ATAN(temp[*,1],temp[*,0]) * !radeg
PLOTS,xy[1,*],xy[0,*],NOCLIP=0,/DATA;,COLOR=128
ENDIF
ENDFOR
FREE_LUN,lun
;!p.font=1
;DEVICE, SET_FONT='Helvetica Bold Italic', /TT_FONT
MAP_GRID,lats=[50,55,60,65,70,75,80,85],lons=[-35]
;Contwoyto Lake 65.75 248.75 72.64 298.10
;Fort Simpson 61.76 238.77 67.27 290.66
;Hopen Island 76.51 25.01 71.92 131.66
;Soroya 70.54 22.22 67.30 120.50
;USERSYMBOL,'star',/FILL,SIZE=2
PLOTS,298.10,72.64,PSYM=5 ;cont
XYOUTS,298.10,72.64,'Contwoto Lake',CHARSIZE=2
PLOTS,131.66,71.92,PSYM=5 ;hop
XYOUTS,131.66,71.92,'Hopen Island',CHARSIZE=2
USERSYMBOL,'triangle',/FILL,SIZE=2
PLOTS,290.66,67.27,PSYM=5 ;sim
XYOUTS,290.66,67.27,'Fort Simpson',CHARSIZE=2
PLOTS,120.50,67.30,PSYM=5 ;sor
XYOUTS,120.50,67.30,'Soroya',CHARSIZE=2
;MAP_GRID,lats=[0,90],lons=[298.1,290.66,131.66,120.5],GLINESTYLE=0
;!p.font=0
END
;+
; NAME: GEOSPACE_CONVERT
;
; PURPOSE: Conversion between various geocentric Cartesian
; space physics coordinate systems.
;
; CATEGORY: Coordinate systems
;
; CALLING SEQUENCE: a= GEOSPACE_CONVERT( Position, From, Into )
;
; INPUTS:
; Position a 3 element vector or 3xN element array
; containing the position(s) to be converted.
; Usually in Cartesian coordinates, use
; /FROM_SPHERICAL if in spherical coordinates.
;
; From the names of the initial and final coordinate
; Into systems. Must be a 2-3 letter string from the set
; of 6 allowable coordinate systems:
;
; GEI Geocentric equatorial inertial
; GEO Geographic (rotating earth)
; GSE Geocentric solar ecliptic
; GSM Geocentric solar magnetospheric
; SM Solar magnetic
; MAG Geomagnetic
;
; KEYWORD PARAMETERS: Many of the coordinate systems require time
; information, which can be provided via the
; following keywords. Note that JD starts at
; noon, while UT starts at midnight. Not my idea.
;
; JULIAN_DAY Result = JULDAY(Month, Day, Year)
; UNIVERSAL_TIME
;
; CDF_EPOCH
;
; OUTPUTS: a 3 element vector or 3 x n element array the
; same size as input Position. Usually Cartesian,
; unless the /INTO_SPHERICAL keyword is set, in
; which case it will be (Radius,Latitude,Longitude)
;
; RESTRICTIONS:
;
; PROCEDURE: Construct the appropriate rotation matrices, combine
them
; and apply to the initial coordinates. For a detailed
; description, see:
;
; Space Physics Coordinate Transformations: A User Guide
; M. A. Hapgood, Planetary and Space Science,
; Volume 40, Number 5, pages 711-717, 1992
;
;
; MODIFICATION HISTORY:
; September 8 1993 Written by Brian Jackel, University of Western
Ontario
;
;-
;-----------------------------------------------------------------------------
function GEOSPACE_ROTATION,angle,axis
;
;This is a short utility routine for creating rotation matrices.
;
;zeta is rotation angle in degrees
;axis is rotation axis, a string of 'X','Y' or 'Z'
;
;The basic idea is to
; start with an empty 3x3 matrix
; put a 1 in the diagonal element corresponding to the rotation axis
; put cos(zeta) in the other diagonal elements
; determine the two off-diagonal terms in the same columns and rows as
the cos(zeta) values
; put -sin(zeta) in the term above the diagonal
; put sin(zeta) in the term below
;
;
;element #'s example for a 'Z' rotation
; 0 1 2 cos(zeta) -sin(zeta) 0
; 3 4 5 sin(zeta) cos(zeta) 0
; 6 7 8 0 0 1
;
sin_zeta= sin(angle*!dtor)
cos_zeta= cos(angle*!dtor)
fill= [1.0d0,cos_zeta,cos_zeta,-sin_zeta,sin_zeta]
rmatrix= fltarr(3,3)
CASE STRUPCASE(axis) OF
'X': rmatrix([0,4,8,7,5])= fill
'Y': rmatrix([4,0,8,6,2])= fill
'Z': rmatrix([8,0,4,3,1])= fill
ELSE: MESSAGE,'AXIS parameter must be string of X, Y, or Z'
ENDCASE
RETURN,rmatrix
END
;-------------------------------------------------------------------
function GEOSPACE_CONVERT,Position,From,Into,HELP=help, $
ROTATION_MATRIX=rotation_matrix, $
CDF_EPOCH_DATE=cdf_epoch_date, $
YMDHMS_DATE=ymdhms_date, $
; JULIAN_DATE=julian_date,
UNIVERSAL_TIME=universal_time, $
FROM_SPHERICAL=from_spherical,
INTO_SPHERICAL=into_spherical
;ON_ERROR,2
IF KEYWORD_SET(HELP) THEN BEGIN
DOC_LIBRARY,'Geospace_Convert'
return, -1.0
END
IF (N_PARAMS() EQ 0) THEN MESSAGE,"Error- input value Position not
defined"
siz= SIZE(position)
; IF (siz(0) EQ 0) THEN MESSAGE,"Error- input value Position not
defined"
; IF (siz(0) GT 2) THEN MESSAGE,"Error- Position must be a 3 element
vector, or a 3xn element array"
; IF (siz(1) NE 3) THEN MESSAGE,"Error- Position must be a 3 element
vector, or a 3xn element array"
_from= STRMID( STRCOMPRESS(STRUPCASE(from),/REMOVE_ALL) ,0,3 )
_into= STRMID( STRCOMPRESS(STRUPCASE(into),/REMOVE_ALL) ,0,3 )
IF (_from EQ _into) THEN RETURN,position ;indentity transformation
;
;
;
seconds_in_day= 24.0d0 * 60.0d0 * 60.0d0
seconds_in_century= 36525.0d0 * seconds_in_day
CDF_EPOCH,epoch2000,2000,01,01,12,/COMPUTE
CDF_EPOCH,epochunix,1970,01,01,/COMPUTE
CDF_EPOCH,epochmjd,1858,11,17,/COMPUTE
IF KEYWORD_SET(CDF_EPOCH_DATE) THEN BEGIN
CDF_EPOCH,cdf_epoch_date,year,month,day,hour,minute,second,/BREAKDOWN
IF (year LT 1970) OR (year GT 2100) THEN MESSAGE,'Note- year
'+STRING(year)+' is outside expected range',/INFORMATIONAL
epoch= cdf_epoch_date
ENDIF ELSE IF KEYWORD_SET(YMDHMS_DATE) THEN BEGIN
ymdhms= ymdhms_date
CASE N_ELEMENTS(ymdhms) OF
1:CDF_EPOCH,epoch,ymdhms[0],/COMPUTE
2:CDF_EPOCH,epoch,ymdhms[0],ymdhms[1],/COMPUTE
3:CDF_EPOCH,epoch,ymdhms[0],ymdhms[1],ymdhms[2],/COMPUTE
4:CDF_EPOCH,epoch,ymdhms[0],ymdhms[1],ymdhms[2],ymdhms[3],/COMPUTE
5:CDF_EPOCH,epoch,ymdhms[0],ymdhms[1],ymdhms[2],ymdhms[3],ymdhms[4],/COMPUTE
6:CDF_EPOCH,epoch,ymdhms[0],ymdhms[1],ymdhms[2],ymdhms[3],ymdhms[4],ymdhms[5],/COMPUTE
ELSE:MESSAGE,'Error- keyword parameter YMDHMS must have between 1
(just year) and 6 (year, month, day, hour, minute, second) elements'
ENDCASE
ENDIF ELSE BEGIN
MESSAGE,'Warning- no time provided, using current system
time',/INFORMATIONAL
st= SYSTIME(1)*1.0d3 ;convert from second to milliseconds
epoch= epochunix + st
ENDELSE
; Determine the modified julian date (since 00:00 UT 17 November 1858)
; and time of day in hours. Then calculate "T_0", the time in julian
; centuries from 12:00 UT on 1 January 2000 (epoch 2000) to the previous
; midnight.
;
mjdate= (epoch - epochmjd)/1d3 ;modified julian date
in seconds
utime= (mjdate MOD seconds_in_day) / 3600.0d0 ;time of day in hours
mjdate= LONG(mjdate / seconds_in_day) ;modified julian date
in days, truncated to integer
T_0= (mjdate - 51544.5d0) / 36525.0d0
; If the keyword is set, convert from spherical to
; Cartesian coordinates. Otherwise just make a
; working copy of the position vector.
;
new_position= position
IF KEYWORD_SET(FROM_SPHERICAL) THEN BEGIN
R= position(0,*)
lat_rad= position(1,*) *!dtor
lon_rad= position(2,*) *!dtor
new_position(0,0)= R * cos(lat_rad) * cos(lon_rad)
new_position(1,0)= R * cos(lat_rad) * sin(lon_rad)
new_position(2,0)= R * sin(lat_rad)
ENDIF
; namelist= ['SM','GEO','GEI','GSE','GSM','SM']
; tlist= [-5,-1,2,3,4]
; w= WHERE
;==================================================================================================
;Determine what transformation ("T") matrices are required. This is
just a 6x6
;lookup table, from which we get a sequence of transformations. The
number indicates
;the transformation matrix (using Hapgood's convention), with 0 just a
place keeper.
;
CASE _from OF
'GEI': matrices= [ [0,0,0], [1,0,0], [2,0,0],
[3,2,0], [4,3,2], [5,1,0] ]
'GEO': matrices= [ [-1,0,0,0], [0,0,0,0], [2,-1,0,0],
[3,2,-1,0], [4,3,2,-1], [5,0,0,0] ]
'GSE': matrices= [ [-2,0,0], [1,-2,0], [0,0,0],
[3,0,0], [4,3,0], [5,1,-2] ]
'GSM': matrices= [ [-2,-3,0,0], [1,-2,-3,0], [-3,0,0,0],
[0,0,0,0], [4,0,0,0], [5,1,-2,-3] ]
'SM': matrices= [ [-2,-3,-4,0,0], [1,-2,-3,-4,0], [-3,-4,0,0,0],
[-4,0,0,0,0], [0,0,0,0,0], [5,1,-2,-3,-4] ]
'MAG': matrices= [ [-1,-5,0,0,0], [-5,0,0,0,0], [2,-1,-5,0,0],
[3,2,-1,-5,0], [4,3,2,-1,-5], [0,0,0,0,0] ]
ELSE:MESSAGE,'Parameter From must be one of GEI, GEO, GSE, GSM, SM,
MAG'
ENDCASE
CASE _into OF
'GEI': matrices= REFORM( matrices(*,0) )
'GEO': matrices= REFORM( matrices(*,1) )
'GSE': matrices= REFORM( matrices(*,2) )
'GSM': matrices= REFORM( matrices(*,3) )
'SM' : matrices= REFORM( matrices(*,4) )
'MAG': matrices= REFORM( matrices(*,5) )
ELSE:MESSAGE,'Parameter Into must be one of GEI, GEO, GSE, GSM, SM,
MAG'
ENDCASE
matrices= matrices( WHERE(matrices NE 0) ) ;throw away any
zeros
;--------------------------------------------------------------------------------------------
;============================================================================================
;This is where any required rotation matrices are actually constructed.
;Make a stack of 6 3x3 matrices, but only calculate the necessary ones.
;Note that some matrices may require partial or full evaluation of
;other matrices, so group some calculations together depending on what
;is required. The result is a little cryptic, but reasonably efficient.
;
;Make a stack for 6 3x3 rotation matrices, only the top five of which
will be used.
;Leave the 0th one empty to keep consistent with notation in Hapgood
tstack= FLTARR(3,3,6)
w1= TOTAL( ABS(matrices) EQ 1 )
w2= TOTAL( ABS(matrices) EQ 2 )
w3= TOTAL( ABS(matrices) EQ 3 )
w4= TOTAL( ABS(matrices) EQ 4 )
w5= TOTAL( ABS(matrices) EQ 5 )
IF (w1 OR w3 or w4) THEN BEGIN
theta= 100.461d0 + 36000.770d0*T_0 + 15.04107d0*utime
;Greenwich mean sidereal time
tstack(0,0,1) = GEOSPACE_ROTATION(theta,'Z') ;"T1" in Hapgood,
equation 2
ENDIF
IF (w2 OR w3 OR w4) THEN BEGIN
ecliptic_obliquity= 23.439d0 - 0.013d0*T_0
mean_anomaly= 357.528d0 + 35999.050d0*T_0 + 0.04107d0*utime
mean_longitude= 280.460d0 + 36000.772d0*T_0 + 0.04107d0*utime
ecliptic_longitude= mean_longitude + (1.915d0 -
0.0048d0*T_0)*sin(mean_anomaly*!dtor) + 0.020*sin(2*mean_anomaly*!dtor)
Ez= GEOSPACE_ROTATION(ecliptic_longitude,'Z')
Ex= GEOSPACE_ROTATION(ecliptic_obliquity,'X')
tstack(0,0,2) = Ez##Ex
ENDIF
IF (w3 or w4 or w5) THEN BEGIN
dyear= (mjdate-46066.0d0)/365.25 ;should be centered on 1985
phi= 78.7646 + 0.0351936*dyear -0.000128034*dyear^2
lambda= 289.097 - 0.0363523*dyear -0.00142802*dyear^2
IF (w5) THEN BEGIN
Ez= GEOSPACE_ROTATION(lambda,'Z')
Ey= GEOSPACE_ROTATION(phi-90.0,'Y')
tstack(0,0,5) = Ey##Ez
ENDIF
ENDIF
IF (w3 or w4) THEN BEGIN
phi_rad= phi*!dtor
cos_phi_rad= cos(phi_rad)
lambda_rad= lambda*!dtor
Qg= [ cos_phi_rad*cos(lambda_rad), cos_phi_rad*sin(lambda_rad),
sin(phi_rad) ]
Qe= REFORM(tstack[*,*,2]) ## ( TRANSPOSE(REFORM(tstack[*,*,1]))##Qg
) ;equation (7)
IF (w3) THEN BEGIN
psi= ATAN( Qe(1),Qe(2) )*!radeg
tstack(0,0,3) = GEOSPACE_ROTATION(-psi,'X')
ENDIF
IF (w4) THEN BEGIN
mu= ATAN( Qe(0), SQRT( Qe(1)^2+Qe(2)^2 ) )*!radeg
tstack(0,0,4)= GEOSPACE_ROTATION(-mu,'Y')
ENDIF
ENDIF
;
;Combine the "T" matrices, do inversions (transpositions) on
;those matrices that need it
;
nstack= N_ELEMENTS(matrices)
T= fltarr(3,3)
T([0,4,8])= [1.0,1.0,1.0] ;set diagonal elements equal to 1
(identity matrix)
FOR i=nstack-1,0,-1 DO BEGIN
temp= REFORM( tstack(*,*,ABS(matrices(i))) )
IF ( matrices(i) LT 0 ) THEN temp= TRANSPOSE(temp)
T= temp##T
ENDFOR
;
;Apply rotation
;
new_position= T ## new_position
rotation_matrix= T
;
;Maybe convert back from Cartesian into spherical
;
IF KEYWORD_SET(INTO_SPHERICAL) THEN BEGIN
X= new_position(*,0)
Y= new_position(*,1)
Z= new_position(*,2)
l= SQRT( X^2 + Y^2 )
new_position(0,0)= SQRT( l^2 + Z^2 )
new_position(0,1)= ATAN(z,l) * !radeg
new_position(0,2)= ATAN(y,x) * !radeg ;Hapgood uses an ArcCos, but
this should work
ENDIF
RETURN,new_position
END
PRO GEOSPACE_CONVERT_TEST
;Formats and units:
; Day/Time format: YY/MM/DD HH.HHHHH
; Degrees/Hemisphere format: Decimal degrees with 2 place(s).
; Longitude 0 to 360, latitude -90 to 90.
; Distance format: Earth radii with 5 place(s).
;
;polar
; Time GEI (RE) GEO
(RE) GM (RE) GSE
(RE) GSM (RE) SM (RE)
;yy/mm/dd hh.hhhhh X Y Z X
Y Z X Y Z X
Y Z X Y Z X
Y Z
;
;00/01/01 0.00000 1.89412 1.02594 6.89192 0.68262
-2.04311 6.89192 0.85108 0.00458 7.17039 -3.30372
2.49712 5.91512 -3.30372 0.82850 6.36693 -0.19481
0.82850 7.17039
;00/01/01 4.95000 1.85155 0.17499 -2.04050 -1.82577
-0.35419 -2.04050 0.14029 -1.84442 -2.04962 0.96501
1.70908 -1.94173 0.96501 1.83839 -1.81978 -0.20466
1.83839 -2.04962
;00/01/01 9.90000 7.25976 1.80121 3.82529 -4.29675
6.12262 3.82529 -7.74501 -2.15036 2.44379 -1.82681
7.70989 2.79316 -1.82681 7.99081 1.84148 -0.86959
7.99081 2.44379
;00/01/01 14.85000 5.14900 1.70083 7.19941 3.11398
4.43939 7.19941 -4.50089 4.35341 6.48277 -3.41217
5.86883 5.92878 -3.41217 5.95795 5.83921 -1.92691
5.95795 6.48277
;00/01/01 19.80000 -0.41009 0.40680 5.17012 -0.07490
0.57276 5.17012 -1.50754 0.10925 4.97787 -2.46367
0.04846 4.58167 -2.46367 -0.92680 4.48722 -1.19400
-0.92680 4.97787
;00/01/01 0.00000 1.89412 1.02594 6.89192 0.68262
-2.04311 6.89192 0.85108 0.00458 7.17039 -3.30372
2.49712 5.91512 -3.30372 0.82850 6.36693 -0.19481
0.82850 7.17039
CDF_EPOCH,epoch,2000,01,01,00,00,00,/COMPUTE
xgei= [1.89412 , 1.02594 , 6.89192 ]
xgeo= [0.68262 , -2.04311 , 6.89192 ]
xmag= [0.85108 , 0.00458 , 7.17039 ]
xgse= [-3.30372 , 2.49712 , 5.91512 ]
xgsm= [-3.30372 , 0.82850 , 6.36693 ]
xsm= [-0.19481 , 0.82850 , 7.17039 ]
PRINT,'MAG to GEO
',SQRT(TOTAL((xgeo-GEOSPACE_CONVERT(xmag,'mag','geo',CDF_EPOCH=epoch))^2))
PRINT,'GEO to GEI
',SQRT(TOTAL((xgei-GEOSPACE_CONVERT(xgeo,'geo','gei',CDF_EPOCH=epoch))^2))
PRINT,'GEI to GSE
',SQRT(TOTAL((xgse-GEOSPACE_CONVERT(xgei,'gei','gse',CDF_EPOCH=epoch))^2))
PRINT,'GSE to GSM
',SQRT(TOTAL((xgsm-GEOSPACE_CONVERT(xgse,'gse','gsm',CDF_EPOCH=epoch))^2))
PRINT,'GSM to SM
',SQRT(TOTAL((xsm-GEOSPACE_CONVERT(xgsm,'gsm','sm',CDF_EPOCH=epoch))^2))
PRINT,'SM to MAG
',SQRT(TOTAL((xmag-GEOSPACE_CONVERT(xsm,'sm','mag',CDF_EPOCH=epoch))^2))
PRINT,'GSM to GEO
',SQRT(TOTAL((xgeo-GEOSPACE_CONVERT(xgsm,'gsm','geo',CDF_EPOCH=epoch))^2))
PRINT,' '
;00/01/01 4.95000 1.85155 0.17499 -2.04050 -1.82577
-0.35419 -2.04050 0.14029 -1.84442 -2.04962 0.96501
1.70908 -1.94173 0.96501 1.83839 -1.81978 -0.20466
1.83839 -2.04962
CDF_EPOCH,epoch,2000,01,01,04,57,00,/COMPUTE
xgei= [ 1.85155 , 0.17499 , -2.04050]
xgeo= [-1.82577 , -0.35419 , -2.04050]
xmag= [0.14029 , -1.84442 , -2.04962 ]
xgse= [ 0.96501 , 1.70908 , -1.94173]
xgsm= [ 0.96501 , 1.83839 , -1.81978]
xsm= [-0.20466 , 1.83839 , -2.04962 ]
PRINT,'MAG to GEO
',SQRT(TOTAL((xgeo-GEOSPACE_CONVERT(xmag,'mag','geo',CDF_EPOCH=epoch))^2))
PRINT,'GEO to GEI
',SQRT(TOTAL((xgei-GEOSPACE_CONVERT(xgeo,'geo','gei',CDF_EPOCH=epoch))^2))
PRINT,'GEI to GSE
',SQRT(TOTAL((xgse-GEOSPACE_CONVERT(xgei,'gei','gse',CDF_EPOCH=epoch))^2))
PRINT,'GSE to GSM
',SQRT(TOTAL((xgsm-GEOSPACE_CONVERT(xgse,'gse','gsm',CDF_EPOCH=epoch))^2))
PRINT,'GSM to SM
',SQRT(TOTAL((xsm-GEOSPACE_CONVERT(xgsm,'gsm','sm',CDF_EPOCH=epoch))^2))
PRINT,'SM to MAG
',SQRT(TOTAL((xmag-GEOSPACE_CONVERT(xsm,'sm','mag',CDF_EPOCH=epoch))^2))
PRINT,'GSM to GEO
',SQRT(TOTAL((xgeo-GEOSPACE_CONVERT(xgsm,'gsm','geo',CDF_EPOCH=epoch))^2))
RETURN
END