小编Rob*_*ans的帖子

我正在尝试将具有高分辨率(25米)和分类数据(1到13)的森林覆盖栅格重新采样为RasterLayer具有较低分辨率(~1 km)的新光栅.我的想法是将森林覆盖数据与其他低分辨率栅格数据结合起来:

1. 我试过raster::resample(),但由于数据是绝对的,我丢失了很多信息:

   summary(as.factor(df$loss_year_mosaic_30m))
     0       1   2   3  4   5   6   7  8   9   10  11   12  13
   3777691  65  101 50 151 145 159 295 291 134 102 126 104  91 
   

   Run Code Online (Sandbox Code Playgroud)

   如您所见,新栅格具有所需的分辨率,但也有很多零.我认为这是正常的,因为我使用了'ngb'选项resample.

2. 第二种策略是使用,raster::aggregate()但我发现难以定义因子整数,因为分辨率的变化不是直截了当的(如分辨率的两倍或相似).

   我的高分辨率栅格具有以下分辨率,我希望它将其聚合0.008333333, 0.008333333 (x, y)到相同程度的分辨率.

   loss_year
   class       : RasterLayer 
   dimensions  : 70503, 59566, 4199581698  (nrow, ncol, ncell)
   resolution  : 0.00025, 0.00025  (x, y)
   extent      : -81.73875, -66.84725, -4.2285, 13.39725  (xmin, xmax, ymin, ymax)
   coord. ref. : +proj=longlat …

   Run Code Online (Sandbox Code Playgroud)

以每月的温度数据多年的光栅文件,它具有通过连接访问的名称names(object)按以下格式"Jan.1981",两年与下面的代码工作"Feb.1981"等(例如文件位置 -添加的所有文件太大了

使用以下代码读入并写入NetCDF:

#Load Packages
library(raster)
library(ncdf4)

#Read in temperature files
r1 <- brick('TavgM_1981.grd')
r2 <- brick('TavgM_1982.grd')

#stack them together 
TempStack = stack(r1, r2)

#set the coordinate system (as it was missing)
crs(TempStack) <- ('+proj=lcc +lat_1=53.5 +lat_2=53.5 +lat_0=46.834 +lon_0=5 +x_0=1488375 +y_0=-203375 +datum=WGS84 +to_meter=2500 +no_defs +ellps=WGS84 +towgs84=0,0,0')

#reproject to get in lat/lon instead of meters
TempStack<-projectRaster(TempStack, crs=CRS("+init=epsg:4326"))

#Extract monthly data names to assign to netCDf later
names <- names(TempStack)

#write the raster file to NetCDF
writeRaster(TempStack, "Temp.nc", overwrite=TRUE, …

我尝试使用 getData （来自栅格包）读取 DEM 栅格，然后将 RasterLayer 转换为 SpatRaster （terra 包）。第一步成功了，但第二步失败了。

library(raster)
library(terra)
 
(alt <- getData('alt', country='DEU', mask=T))
#class      : RasterLayer 
#dimensions : 960, 1116, 1071360  (nrow, ncol, ncell)
#resolution : 0.008333333, 0.008333333  (x, y)
#extent     : 5.8, 15.1, 47.2, 55.2  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : D:/dummy/DEU_msk_alt.grd 
#names      : DEU_msk_alt 
#values     : -179, 2701  (min, max)

altT <- rast(alt)            
# rast is supposed to be able to read RasterLayer, but it fails. Why?
# Error …

我正在使用 \xe2\x80\x9cterra\xe2\x80\x9d 包处理 SpatRaster 列表。我已经从脚本创建了栅格列表，并保存了 R 环境。但是，当我从另一个脚本加载 SpatRasters 列表时，遇到以下消息错误：

Error: external pointer is not valid\n

这是一个可重现的示例：

library(terra)\nx <- terra::rast(xmin=-110, xmax=-80, ymin=40, ymax=70, ncols=30, nrows=30)\nvalues(x) <- 1:ncell(x)\nr <- c(x, x, x, x)\nr <- list(r, r, r, r)\nsave(r, file = "test.Rdata")\n\nrm(list=ls(all=TRUE))\n\nload("test.Rdata")\nr\n#[[1]]\n#class       : SpatRaster \n#Error: external pointer is not valid\n

您能否提供解决此问题的指导？\xc2\xa0任何帮助将不胜感激。

我遇到了一个不寻常的结果，R 中的栅格没有数据值。下面的代码 - 我有一个没有导入的数据值的栅格 (-9999)。QGIS 也读取 nodatavalue，arcgis。R 在读入 geotiff 时，将 nodata 值分配给 -INF。

我不知道为什么。但是我决定尝试从头开始制作一个 - 并且结果相同。我的过程有什么问题？如何确保 R 正确读取 nodatavalues？

注意：我下面的示例是创建的栅格/ geotiff。但我实际上是在导入一个组织生产的相当大的 geotiff。所以我无法控制它们的编写方式，但如果需要，我可以要求他们调整标签。

library(raster)
#create a raster from the matrix
myRaster1 <- raster(nrow=4, ncol=4)

#assign some random data to the raster
myRaster1[]<- 1:ncell(myRaster1)

myRaster1[5] <- -9999

#ensure the data have some decimals
myRaster1[2] <- 34.5

#assign no data value to raster
myRaster1@file@nodatavalue <- -9999

#make sure it worked
NAvalue(myRaster1)
myRaster1@file@nodatavalue

#view attributes of the raster
myRaster1


#write out raster
#write the …

我在栅格化 shapefile 以在 0.5*0.5 网格上生成点时遇到问题。shapefile 表示全球珊瑚礁对综合威胁的风险级别（低 0、中 100、高 1000、极高 1500）的分类。

我从另一个工作正常的示例中提取了代码，但是当我为我的数据尝试它时，我从 plot 函数中什么也得不到。有关 shapefile 和我的代码的链接，请参见下文：

珊瑚礁面临风险：全球综合威胁

# Read shapefile into R
library(rgdal)
library(raster)    

int.threat.2030 <- readOGR(dsn = "Global_Threats/Integrated_Future", 
                           layer = "rf_int_2030_poly")

## Set up a raster "template" for a 0.5 degree grid
ext <- extent(-110, -50, 0, 35)
gridsize <- 0.5
r <- raster(ext, res=gridsize)

## Rasterize the shapefile
rr <- rasterize(int.threat.2030, r)

## Plot raster
plot(rr)

我可能会出错的任何想法？是shapefile本身的问题吗？

请和谢谢！

我有一个函数,它以光栅砖对象的形式读取多波段图像,迭代通过各种计算的波段,然后将光栅写为新的.tif.所有这一切都很好,但新图像文件的文件大小大约是其四倍(我假设因为原始图像有4个波段).我想知道在writeRaster()函数中是否有一个我不知道的参数,或者如果有其他方法我可以确保输出图像与输入的文件大小基本相同.

原始文件大小为134 MB; 输出范围从471到530 MB左右,具体取决于格式.

简化代码:

library(rgdal)
library(raster)

path = "/Volumes/ENVI Standard Files/"
img = "qb_tile.img"

imageCorrection = function(path, img){
  raster = brick(paste0(path, img))   
  raster = reclassify(raster, cbind(0, NA))  

  for(i in 1:nlayers(raster)){   
    raster[[i]] = raster[[i]] - minValue(raster[[i]]) 
  }
  writeRaster(raster, paste0(path,img,"_process.tif"), format = "GTiff", overwrite=TRUE)
}

我注意到，如果我使用 terra 包从存储在本地硬盘上的文件加载栅格，然后稍后加载工作区，则 SpatRaster 对象没有任何与其关联的数据。有没有办法在保存和加载工作空间时保留与 SpatRaster 对象关联的所有信息？

这是一个示例代码来说明这个问题：

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)

#This produces the following output
r
#class       : SpatRaster 
#dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
#resolution  : 0.008333333, 0.008333333  (x, y)
#extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, 
               ymin, #ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#source      : elev.tif 
#name        : elevation 
#min value   :       141 
#max value   :       547

sources(r)#this works

save.image("delete_if_found.RData")

rm(list = ls())

load("delete_if_found.RData")
r
#which returns …

我正在尝试将与蝙蝠位置 (SpatialPointsDataFrame) 相关的数据叠加到科罗拉多州 (SpatialPolygonsDataFrame) 上。两个对象的CRS不同：

crs(colorado)
#CRS arguments:
# +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
crs(BatRange_data)
#CRS arguments:
# +proj=utm +zone=13 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 

当我重新投影蝙蝠数据时，CRS确实发生了变化，但程度不受影响。

前：

BatRange_data
#class       : SpatialPointsDataFrame 
#features    : 2456 
#extent      : 139996.3, 748812, 4094998, 4535103  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=utm +zone=13 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
#variables   : 17

后：

geo_proj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
crs(BatRange_data) = geo_proj

BatRange_trnsfrmd = spTransform(BatRange_data,geo_proj)
BatRange_trnsfrmd
#class       : SpatialPointsDataFrame 
#features    : 2456 
#extent      : …

这基本上是在这个问题https://gis.stackexchange.com/questions/95481/in-r-set-na-cells-in-one-raster-where-another-raster-has-values和此处为大栅格图层中的 R-Chaging 特定像元值。第一个问题的答案没有解决我的问题，因为我没有另一个可以与overlay或 一起使用的对象calc。重新分类是唯一有效的方法。

就我而言，我有一个 300 MB 的光栅文件。我正在应用一个简单的操作，只是尝试用 NA 替换栅格中等于某个数字的所有值。我有大约 4 GB 的可用 RAM，但似乎我无法完成此操作，因为我收到错误消息“无法分配大小为 4.6 GB 的向量”。我什至尝试将我的内存大小设置为 16 GB，但是我得到了同样的错误，只是说无法分配大小为 9.2 GB 的向量。我尝试了以下两个选项：

r[r==5]=NA 
values(r)[values(r)==5]

奇怪的是，即使是像 table(values(r)) 这样的简单操作也会出现同样的错误，例如 ArcMap 可以在几秒钟内创建这个表。我已经解决了我的问题，但我想知道为什么内存使用效率如此低下，以及如何防止或避免它？为什么raster需要多达 9 GB 来处理 300 MB 的文件？这是这个包的限制还是 R 的限制？