星图中的球极平面投影

将球面投影到平面上,有多种投影方式。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")
}