从栅格中提取值时返回的 NA答案

【问题标题】：NAs returned when extracting values from a raster从栅格中提取值时返回的 NA
【发布时间】：2018-12-01 02:35:36
【问题描述】：

我正在尝试使用空间点文件从显示熊密度的栅格中提取值，我的脚本如下：

## set the working directory                                                         
##

setwd("~/Masters/Research Project/R Scripts/Focal_Raster")

## install packages to use                                                           
##

library(raster)
library (rgdal)

## import the raster and the csv file that will be used in the habitat 
## location      ##
## extraction                                                                       

r <- raster("bear_12zero2.tif")
r
plot(r, main = "Bear Density")

hist(r[])
head(r[],200)

locations<-read.csv("PA_AllBears_1.csv")
plot(locations$X,locations$Y, main="BearID")

## first need to assign a coordinate reference system to the csv file of 
## data        ##
## locations that we want to extract                                                 

crs(locations)
library(rgdal)
library(sp)
crs(r)

locationsC = SpatialPoints(cbind(locations$X, locations$Y), 
                       proj4string=CRS("+init=epsg:4326"))
cord.UTM <- spTransform(locationsC, crs(r))
cord.UTM


## Now we just need to extract the raster values from these locations                
##
TRI<-extract(r,cord.UTM,method='simple')

TRI

## need to write CSV file to join to the master datasheet with all varibles 
## in       ##

write.csv(TRI,'beardensity.csv')

当我运行脚本时，光栅显示：

> r
class       : RasterLayer 
dimensions  : 2000, 1394, 2788000  (nrow, ncol, ncell)
resolution  : 2248.04, 2248.04  (x, y)
extent      : -27831.02, 3105937, 4846432, 9342512  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
+towgs84=0,0,0 
data source : C:\Users\Kate\Documents\Masters\Research Project\R 
Scripts\Focal_Raster\bear_12zero2.tif 
names       : bear_12zero2

但是当我尝试使用我已导入并转换为相同投影的位置 csv 提取值时，我最终得到了所有 NA。我使用的另一个栅格有类似的错误，它具有相同的投影并使用相同的提取点。

我的结果如下：

> setwd("~/Masters/Research Project/R Scripts/Focal_Raster")
> library(raster)
> library (rgdal)
> r <- raster("bear_12zero2.tif")
> r
class       : RasterLayer 
dimensions  : 2000, 1394, 2788000  (nrow, ncol, ncell)
resolution  : 2248.04, 2248.04  (x, y)
extent      : -27831.02, 3105937, 4846432, 9342512  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
+towgs84=0,0,0 
data source : C:\Users\Kate\Documents\Masters\Research Project\R 
Scripts\Focal_Raster\bear_12zero2.tif 
names       : bear_12zero2 

> plot(r, main = "Bear Density")
> locations<-read.csv("PA_AllBears_1.csv")
> plot(locations$X,locations$Y, main="BearID")
> summary(locations)
     BearID             X               Y              P_A     
 Min.   : 2.000   Min.   :18.37   Min.   :65.61   Min.   :0.0  
 1st Qu.: 5.000   1st Qu.:20.11   1st Qu.:66.34   1st Qu.:0.0  
 Median : 8.000   Median :21.09   Median :66.53   Median :0.5  
 Mean   : 8.273   Mean   :20.78   Mean   :66.51   Mean   :0.5  
 3rd Qu.:12.000   3rd Qu.:21.46   3rd Qu.:66.74   3rd Qu.:1.0  
 Max.   :15.000   Max.   :22.35   Max.   :67.08   Max.   :1.0  
> str(locations)
 'data.frame':  5500 obs. of  4 variables:
 $ BearID: int  2 2 2 2 2 2 2 2 2 2 ...
 $ X     : num  19.5 19 18.8 19.4 19.7 ...
 $ Y     : num  66.4 66.1 66.2 66.3 66.2 ...
 $ P_A   : int  0 0 0 0 0 0 0 0 0 0 ...
> crs(locations)
[1] NA
> library(rgdal)
> library(sp)
> crs(r)
CRS arguments:
 +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
 +towgs84=0,0,0 
> locationsC = SpatialPoints(cbind(locations$X, locations$Y), 
+                            proj4string=CRS(("+init=epsg:4326")))
> cord.UTM <- spTransform(locationsC, crs(r))
> cord.UTM
class       : SpatialPoints 
features    : 5500 
extent      : 650979.3, 825051.9, 7280884, 7457682  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
+towgs84=0,0,0 
> TRI<-extract(r,cord.UTM,method='simple')
> TRI
   [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
  [37] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
  [73] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [109] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [145] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [181] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [217] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [253] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [289] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [325] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [361] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [397] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [433] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [469] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA
 [505] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
NA NA NA NA NA NA NA NA NA NA NA NA NA

 [ reached getOption("max.print") -- omitted 4500 entries ]

任何有用的提示将不胜感激

【问题讨论】：

你为什么要投影你的位置，然后在你的提取中使用未转换的点数据？你的提取表达式应该是extract(r, cord.UTM, method = 'simple')吗？
抱歉，JMT2080AD 正在尝试不同的方法来使脚本正常工作，但忘记将其更改回来。刚刚使用转换后的位置再次运行脚本，但仍然遇到与 NA 相同的问题。
您应该使用更正后的脚本更新您的帖子。
运行plot(r); plot(cord.UTM, add = T) 时看到了什么？点是否绘制在有值区域的栅格上？
当我绘制绘图时，我可以看到我的位置点没有通过值传递到栅格上。

标签： r extract raster

【解决方案1】：

下面的脚本是否给出了预期的结果？如果不查看您拥有的数据，很难知道。我正在使用栅格投影来重新投影您的坐标。

我也关心这个表达式：

CRS("+proj=utm")

您没有定义您想要的 UTM。我不认为有“UTM”，它是“UTM Zone 10 North”，或者类似的东西。有关更多信息，请参阅http://spatialreference.org/。如果你不确定，试试这个网站：http://projfinder.com/你可以把地图放在数据应该在的地方，输入一个坐标对，它会告诉你最适合这些坐标的投影。

## set the working directory                                                         
##

setwd("~/Masters/Research Project/R Scripts/Focal_Raster")

## attach packages to use                                                           
##

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

## import the raster and the csv file that will be used in the habitat 
## location      ##
## extraction                                                                       

r <- raster("bear_12zero2.tif")
plot(r, main = "Bear Density")

hist(r[])
head(r[],200)

locations<-read.csv("locations.csv")
plot(locations$X,locations$Y, main="BearID")

## first need to assign a coordinate reference system to the csv file of
## locations that we want to extract                                                 

locationsC <- SpatialPoints(cbind(locations$X, locations$Y), 
                            proj4string=CRS("+init=epsg:4326"))
cord.UTM <- spTransform(locationsC, crs(r))

## Now we just need to extract the raster values from these locations                
##
TRI <- extract(r, cord.UTM, method='simple')

## need to write CSV file to join to the master datasheet with all varibles 
## in       ##

write.csv(TRI,'beardensity.csv')

【讨论】：

我已经尝试了使用“+init=epsg:32633”而不是“+proj=utm”的建议，因为这是我正在查看的区域的投影，但我仍然得到了 NAs从我提取值时开始
@kate 你能从 csv 中发布你的坐标的最大值和最小值吗？
这些是来自转换后的 CSV > cord.UTM 类：SpatialPoints 特征：5500 范围：650979.3、825051.9、7280884、7457682（xmin、xmax、ymin、ymax）坐标。参考。 : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
来自未转换的csv如下
X 最大值 = 22.38，最小值 = 18.37； Y 最大值= 67.08，最小值= 65.61

【解决方案2】：

我在这个问题上遇到了同样的问题，我非常关注光栅是问题所在。在@JMT2080AD 发表评论后，我绘制并意识到我的观点与我的栅格不一致。所以我改变了设置我的观点的方式：请注意“sites”和“stacked”是我之前导入的积分和rasterStack

coordinates(sites) <- ~Longitude + Latitude

crs(sites) <- crs(stacked)

values <- extract(stacked, sites)

这立即解决了问题，所以也许这也会帮助其他人，以及将来遇到类似问题的人。谢谢。

【讨论】：