如何在R中将阿拉斯加空间点投影到重缩放平移后的美国地图

我看到你已经成功用fixup和fix1函数调整了阿拉斯加和夏威夷的位置，但阿拉斯加的点无法正常显示，核心问题在于点数据的投影处理和变换逻辑与地图不匹配。下面是具体的解决思路和修正后的完整代码：

问题根源分析 #

你的my_points是WGS84经纬度格式（对应epsg:4326），但你直接给它设置了epsg:2163的投影，没有做正确的投影转换；更关键的是，重定位后的阿拉斯加地图经过了旋转、缩放和平移操作，而你的阿拉斯加点没有同步执行这些变换，导致它们还停留在原始的高纬度位置，无法和重定位后的地图区域重合。

修正步骤 #

1. 正确设置点数据的初始投影（经纬度epsg:4326），并转换到epsg:2163投影（和地图调整时使用的投影一致）。
2. 筛选出阿拉斯加的点，对其应用和地图阿拉斯加完全相同的fix1变换（旋转、缩放、平移）。
3. 合并经变换的阿拉斯加点和本土点，再一起绘制到重定位后的地图上。

修正后的完整代码 #

require(maptools)
require(rgdal)
require(sp)

# 保留你原有的fixup和fix1函数
fixup <- function(usa,alaskaFix,hawaiiFix){
 alaska=usa[usa$STATE_NAME=="Alaska",]
 alaska = fix1(alaska,alaskaFix)
 proj4string(alaska) <- proj4string(usa)
 hawaii = usa[usa$STATE_NAME=="Hawaii",]
 hawaii = fix1(hawaii,hawaiiFix)
 proj4string(hawaii) <- proj4string(usa)
 usa = usa[! usa$STATE_NAME %in% c("Alaska","Hawaii"),]
 usa = rbind(usa,alaska,hawaii)
 return(usa)
}

fix1 <- function(object,params){
 r=params[1];scale=params[2];shift=params[3:4]
 object = elide(object,rotate=r)
 size = max(apply(bbox(object),1,diff))/scale
 object = elide(object,scale=size)
 object = elide(object,shift=shift)
 object
}

# 下载并读取美国州界shapefile
setwd(tempdir())
download.file("https://dl.dropbox.com/s/wl0z5rpygtowqbf/states_21basic.zip?dl=1", "usmapdata.zip", method = "curl")
unzip("usmapdata.zip")

us = readOGR(dsn = "states_21basic",layer="states")

# 转换投影并调整阿拉斯加、夏威夷位置
usAEA = spTransform(us,CRS("+init=epsg:2163"))
usfix = fixup(usAEA,c(-35,1.5,-2800000,-2600000),c(-35,1,6800000,-1600000))
plot(usfix, main = "重定位后的美国地图（含所有点位）")

# 处理你的点数据
my_points= structure(c(44.334567, 44.209571, 44.049845, 44.752622, 44.791511, 44.391792, 44.540124, 46.075669, 46.007611, 45.836781, 46.595384, 45.703999, 45.47594, 44.715117, 44.385954, 42.804004, 43.349842, 43.252618, 42.891499, 56.212978, 55.391604, 55.395194, 60.476929, 61.709743, 62.76729, 62.346439, 59.93483, 64.472359, 64.88569, -122.047007, -122.256733, -123.42621, -122.128965, -122.578974, -122.497582, -122.435915, -121.998698, -122.347319, -122.466208, -122.459556, -123.755405, -123.725121, -123.887335, -123.831778, -123.610907, -122.728942, -123.026172, -124.070652, -131.638358, -131.195583, -132.408627, -151.081668, -149.231938, -149.693379, -150.019201, -158.190006, -146.926265, -147.249648), .Dim = c(29L, 2L))

# 转换为空间点对象，初始投影设置为epsg:4326（经纬度）
my_NEW_S <- as.data.frame(my_points)
colnames(my_NEW_S) <- c("y", "x")
coordinates(my_NEW_S) <- ~x+y 
proj4string(my_NEW_S) <- CRS("+init=epsg:4326") 

# 将点转换到epsg:2163投影（与地图调整时的投影一致）
my_NEW_S_AEA <- spTransform(my_NEW_S, CRS("+init=epsg:2163"))

# 筛选阿拉斯加的点（根据数据特征，纬度>50的为阿拉斯加点）
ak_points <- my_NEW_S_AEA[my_NEW_S_AEA$y > 50, ]
contiguous_points <- my_NEW_S_AEA[my_NEW_S_AEA$y <= 50, ]

# 对阿拉斯加点应用和地图完全相同的fix1变换
ak_points_fixed <- fix1(ak_points, c(-35,1.5,-2800000,-2600000))
proj4string(ak_points_fixed) <- proj4string(my_NEW_S_AEA)

# 合并所有点位
all_points_fixed <- rbind(contiguous_points, ak_points_fixed)

# 在重定位后的地图上绘制所有点
points(all_points_fixed, col="red", pch=16, cex=0.5)

关键细节说明 #

• 先把点从经纬度（epsg:4326）转换到epsg:2163投影，确保和调整后的地图处于同一空间坐标系。
• 用纬度筛选阿拉斯加点是基于你的数据特征（所有阿拉斯加点纬度都高于50），如果你的数据有差异，可以根据实际情况调整筛选条件。
• 对阿拉斯加点使用和地图完全一致的fix1参数，保证它们的旋转、缩放、平移操作和地图同步，这样就能完美匹配重定位后的阿拉斯加区域。

内容的提问来源于stack exchange，提问作者Ali