r - 将空间坐标分配给 R 中的数组

Question

我已从以下链接下载了数据文本文件:http://radon.unibas.ch/files/us_rn_50km.zip

解压缩后，我使用以下代码行来绘制数据: # 加载库 库(字段)

# function to rotate a matrix (and transpose)
rotate <- function(x) t(apply(x, 2, rev))

# read data
data <- as.matrix(read.table("~/Downloads/us_rn_50km.txt", skip=6))
data[data<=0] <- NA

# rotate data
data <- rotate(data)

# plot data
mean.rn <- mean(data, na.rm=T)
image.plot(data, main=paste("Mean Rn emissions =", sprintf("%.3f", mean.rn)) )

这一切看起来都不错，但我希望能够在经纬度网格上绘制数据。我想我需要将此数组转换为 sp 类对象，但我不知道如何。我知道以下内容(来自网站):“用于投影纬度和经度坐标的投影是用于北美地质十年 (DNAG) map 的投影。投影类型是球形横轴墨卡托投影，基本纬度为零度，引用经度为西经 100 度。使用的比例因子为 0.926，无东偏移或北偏移。经纬度基准面为 NAD27，xy 坐标的单位为米。使用的椭球体为 Clarke 1866。 map 分辨率为50x50km”。但不知道如何处理这些数据。我尝试过:

proj4string(data)=CRS("+init=epsg:4267")
data.sp <- SpatialPoints(data, CRS("+proj=longlat+ellps=clrk66+datum=NAD27") )

但是遇到了各种问题(NA)，从根本上来说，我认为数据的格式不正确。我认为 SpatialPoints 函数需要位置数据(二维)和第三个值数组与这些位置相关联(x、y、z 数据 - 我想我的问题是从我的数据中计算出 x 和 y!)

非常感谢任何帮助!

谢谢， 亚历克斯

Answer 1

相关文件是 ASCII 栅格网格。坐标在此格式中是隐式的；标题描述(通常)左下角的位置，以及网格尺寸和分辨率。在此标题部分之后，由空格分隔的值描述了变量在网格中的变化方式，值按行优先顺序给出。如果您感兴趣，请在文本编辑器中打开它。

您可以使用出色的raster包将此类文件导入到R中，如下所示:

download.file('http://radon.unibas.ch/files/us_rn_50km.zip', 
          destfile={f <- tempfile()})
unzip(f, exdir=tempdir())
r <- raster(file.path(tempdir(), 'us_rn_50km.txt'))

您可以立即绘制它，而无需指定投影:

如果您不想将其转换为另一个 CRS，则不一定需要指定当前坐标系。但由于您确实想要将其转换为 WGS84(地理)，因此您需要首先分配 CRS:

proj4string(r) <- '+proj=tmerc +lon_0=-100 +lat_0=0 +k_0=0.926 +datum=NAD27 +ellps=clrk66 +a=6378137 +b=6378137 +units=m +no_defs'

不幸的是，我不完全确定这个 proj4string 是否正确反射(reflect)了 website that provided the data 中给出的信息。 (如果他们确实以标准格式提供定义，那就太好了)。

分配 CRS 后，您可以使用 projectRaster 投影栅格:

r.wgs84 <- projectRaster(r, crs='+init=epsg:4326')

如果您愿意，可以将其写成您选择的光栅格式，例如:

writeRaster(r.wgs84, filename='whatever.tif')