如何在R中对栅格文件应用地理投影？

Question

我正在尝试从60个栅格文件中创建一个大型的（空间？）数据框，以便在随后的多元统计分析（例如PCA）中使用。

我对在R中处理空间数据还不太熟悉，不理解空间坐标数据是否被保留或约束在数据文件中。

我使用以下代码读取栅格文件并创建了一个栈：

#load libraries
library(raster)
library(sp)
library(rgdal)
library(spatial.tools) 

#set wd
setwd("C:/Users/...../data/")

#get raster files from wd 
files <- list.raster.files(path = getwd(), pattern = ".tif$", recursive =FALSE, return_rasters = FALSE, return_bbox = FALSE)

#create rasterStack
mystack <- stack(files$raster_files)

#read metadata / summary of rasterStack
mystack
#output has "NA" for coordinate reference system. 

为什么我导入文件后，R会丢失空间参考信息？当我在ArcGIS中查看栅格文件时，它们确实具有投影信息。

我尝试手动应用以下方式的投影：

#first define the projection using the proj4 syntax 
projection <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"

#then apply the projection to the rasterStack
projectRaster(files, projection, method = bilinear, filename = "STACKprj")

这给我一个错误，即输入数据("mystack")的参照系统为NA，也就是说，如果没有一个起始参照系统，它就无法工作。

#I try instead to apply the projection in this way: 
STACKprj <- projection(mystack)

那行不通。

有什么技巧吗? 谢谢!

- user3251223

---

projection(mystack) <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs" 这段代码有效吗？ - Mike.Gahan

请使用以下代码代替 projection <- CRS("+proj=merc +a=6378137 +b=6378137 ..."); proj4string(mystack) <- projection - Paulo E. Cardoso

Answer 1

你应该尝试类似这样的内容：

library(raster)
library(sp)
library(rgdal)

#work wd
wd <- "C:/Users/...../data"

#get raster files from work wd 
files <- list.raster.files(path = wd, pattern = ".tif$",
  recursive =FALSE, return_rasters = FALSE, return_bbox = FALSE)

#create rasterStack
mystack <- stack(files$raster_files)
proj4string(mystack) <- CRS("+init=epsg:3857") # OSM Mercator projection

- Paulo E. Cardoso

---

非常感谢你，Paulo！你的两个建议都完美地解决了问题！！ - user3251223

不客气！我建议您给最符合您需求（有用的）答案点赞。 - Paulo E. Cardoso