2

我有一个极地立体网格（尺寸为 6667 x 6667 单元格，大小为顶部：3333500，左侧：-3333500，右侧：3333500，底部：-3333500）。该投影的纬度为真实比例南 -71 度，基准 WGS84。网格间距为 1000 m

我想从中制作一个纬度长的网格，但遇到了麻烦。我可能会完全错误地解决这个问题，但这是我迄今为止所拥有的：

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 的输出不正确。有没有人有什么建议？谢谢！

4

2 回答 2

5

这会有帮助吗？

#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)

每股收益 3031

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)

wgs​​84

还有一个更合理的例子

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))

ggmap

编辑

有关 EPSG:3031 WGS 84 / 南极极地立体图的更多信息，请访问NCIDC 网站或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()

在此处输入图像描述

于 2014-05-23T21:48:25.330 回答
0

通过查看此文档页面，似乎+proj=stere需要 a+lat_0和+lon_0参数。使用

stere <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +datum=WGS84 +units=m"

似乎确实进行了转换（至少它具有非零值）。

proj4string(a)<-CRS(stere) 
z<-spTransform(a,CRS("+proj=longlat +datum=WGS84")) 
summary(z)
# Object of class SpatialPointsDataFrame
# Coordinates:
#   min       max
# x -45 135.00000
# y -90 -48.45867
# Is projected: FALSE 
# proj4string :
# [+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]
# Number of points: 6667
于 2014-05-23T21:49:52.167 回答