r - 将极地立体投影转换为R中的经纬度网格

Question

我有一个极地立体网格(尺寸为 6667 x 6667 单元格，范围为顶部:3333500， 左:-3333500，右:3333500，下:-3333500)。该投影的纬度为真实比例尺以南 -71 度，基准面 WGS84。网格间距1000米

我想由此制作一个经纬度网格，但遇到了麻烦。我的做法可能完全错误，但这是我迄今为止所得到的:

library(rgdal)
x<-seq(-3333500,3333500, length.out=6667)
y<-seq(3333500,-3333500,length.out=6667)
a<-data.frame(x,y)
coordinates(a)= ~x + y

stere <- "+proj=stere +lat_ts=-71 +datum=WGS84 +units=m"
#i have also tried:
#stere <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84  "
proj4string(a)<-CRS(stere) 
spTransform(a,CRS("+proj=longlat +datum=WGS84")) 

spTransform 的输出不正确。有没有人有什么建议？谢谢!

Answer 1

这有什么帮助吗？

#Load packages

kpacks <- c("rgdal", 'ggplot2', 'maptools', 'raster')
new.packs <- kpacks[!(kpacks %in% installed.packages()[ ,"Package"])]
if(length(new.packs)) install.packages(new.packs)
lapply(kpacks, require, character.only=T)
remove(kpacks, new.packs)

#Your GRID limits

x<-seq(-3333500,-3333000, length.out=10)
y<-seq(-3333000,-3333500,length.out=10)
xy <- as.data.frame(expand.grid(x,y))
coordinates(xy)= ~Var1 + Var2
plot(xy, axes = T)

proj.pol <- CRS('+init=epsg:3031')
wgs <- CRS('+init=epsg:4326')
proj4string(xy) <- proj.pol
awgs <- spTransform(xy, wgs)
head(awgs)
SpatialPoints:
          Var1      Var2
[1,] -134.9957 -48.46152
[2,] -134.9962 -48.46184
[3,] -134.9967 -48.46216
[4,] -134.9971 -48.46248
[5,] -134.9976 -48.46280
[6,] -134.9981 -48.46311
Coordinate Reference System (CRS) arguments: +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0

plot(awgs)

还有一个更合理的例子

data(wrld_simpl)
maps <- wrld_simpl[wrld_simpl$NAME %in% c("Argentina", "Chile",
                                          "Brazil",  "Antarctica"), ]
mapsdf <- fortify(maps)
x<-seq(-3433000,3433000, length.out=10)
y<-seq(3433000,-3433000,length.out=10)
xy <- as.data.frame(expand.grid(x,y))
coordinates(xy)= ~Var1 + Var2
proj4string(xy) <- proj.pol
awgs <- spTransform(xy, wgs)
#plot(awgs, axes = T)
awgsdf <- as.data.frame(awgs)

ggplot(maps) +
  geom_path(aes(x=long, y= lat, group = group)) +
  geom_point(data = awgsdf, aes(x=x, y=y)) +
  #coord_polar()
  coord_map("ortho", orientation=c(-40, -20, 10))

编辑

有关 EPSG:3031 WGS 84/南极极地立体投影的更多信息，请访问 NCIDC site或在 remotesensing.org

将剪辑网格添加到范围 您可以定义感兴趣的区域来相应地裁剪网格。

x<-seq(-12400000, 12400000, length.out=50)
y<-seq(-12400000, 12400000,length.out=50)
xy <- as.data.frame(expand.grid(x,y))
coordinates(xy)= ~Var1 + Var2
proj4string(xy) <- proj.pol
awgs <- spTransform(xy, wgs)
plot(awgs, axes = T)

# Create a extent object using raster::extent
ext1 <- extent(matrix(c(-60, 60, -86, -40), byrow = T, nrow=2))
awgs1 <- crop(awgs, ext1) # crop spdf to extent
# Plot it
ggplot(maps) +
  geom_polygon(aes(x=long, y= lat, group = group)) +
  geom_point(data = as.data.frame(coordinates(awgs)),
             aes(x=Var1, y=Var2), size = 1, colour = 'grey60') +
  geom_point(data = as.data.frame(coordinates(awgs1)),
             aes(x=Var1, y=Var2)) +
  #coord_polar()
  coord_map("ortho", orientation=c(-60, -20, 10)) +
  theme_bw()