星图中的球极平面投影

将球面投影到平面上，有多种投影方式。R中的mapproj程序包包含了大部分的投影，适用于地图的投影变换。 而星图中最常用的投影方式为球极平面投影（Stereographic projection）， 这种投影方式有时也用于地图中。

这里介绍在R中如何实现球极平面投影。

定义函数stereo_proj

  • lambda： 其中lambda为赤经
  • phi： 为赤纬
  • lambda0 =0： lambda0为地图（星图）的投影中心点 赤经
  • phi0 =0： phio为地图（星图）投影中心点的 赤纬
  • R =1： 为比例半径， 一般取值为1

函数定义如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
stereo_proj <-function(lambda,phi, lambda0 =0, phi0 =0, R =1){
   lambda <- lambda/(180/pi)
   phi    <- phi/(180/pi)
   lambda0 <- lambda0/(180/pi)
   phi0   <- phi0/(180/pi)
   k =2*R/(1+sin(phi0)*sin(phi)+cos(phi0)*cos(phi)*cos(lambda - lambda0))
   x = k*cos(phi)*sin(lambda - lambda0)
   y = k*(cos(phi0)*sin(phi)-sin(phi0)*cos(phi)*cos(lambda - lambda0))   
   return(list(x = x, y = y))
} 

### 公式参见 ### Source
###  
### 下面用地图投影变换举例

 library(maps)
 ###
 aaa <- map()  
 ### 获取大比例尺世界地图中的座标
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
plot(stereo_proj(lambda = aaa$x, phi = aaa$y), 
      type ="l", xlim =c(-6,6), ylim =c(-6,6)) 

### 添加经度线
####Adding longitudinal Lines
long =c(-180, -150, -120, -90, -60, -30, 0, 30, 60, 90, 120, 150)
lat  =-90:90
for(i in long){
   temp <- stereo_proj(lambda = i, phi = lat)
   points(temp$x, temp$y, type ="l")
}
### 添加纬度线
### Adding latitudinal lines### 
plot(0, xlim =c(-2, 2), ylim = c(-2, 2), pch = "")
lat =c(-90, -60, -30, 0, 30, 60, 90)
long  =-150:150

for(i in lat){
   temp <- stereo_proj(lambda = long,phi = i)
   points(temp$x, temp$y, type ="l")
}