从栅格中提取值时返回的NA [英] NAs returned when extracting values from a raster
本文介绍了从栅格中提取值时返回的NA的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我正在尝试使用空间点文件从显示熊的密度的栅格中提取值，我的脚本如下:
I am trying to extract values from a raster showing density of bears using a spatial points file my script is as follows:
## 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')
运行脚本时，栅格显示:
When I run the script the the raster shows:
> 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.我正在使用的另一个栅格具有类似的错误，该栅格具有相同的投影和相同的提取点.
but when I try to extract values using the locations csv I have imported and transformed to the same projection I end up with all NAs. I have a similar error with another raster I am using which has the same projection and using the same extraction points.
我的结果如下:
> 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 ]
任何有用的技巧将不胜感激
Any useful tips would be greatly appreciated
推荐答案
下面的脚本是否给出了预期的结果?不检查您的数据就很难知道.我正在使用栅格的投影来重新投影您的坐标.
Does the script below give the expected results? It's hard to know without reviewing the data you have. I am using the raster's projection to re-project your coordinates.
我也很关心这个表达:
CRS("+proj=utm")
您没有定义所需的UTM.我认为没有"UTM"，"UTM Zone 10 North"之类的东西.有关更多信息，请参见 http://spatialreference.org/.如果不确定，请尝试以下站点: http://projfinder.com/您可以将地图放在哪里输入数据，然后输入一个坐标对，它会告诉您这些坐标最适合的投影是什么.
You are not defining what UTM you want. I don't think there is a "UTM", it's "UTM Zone 10 North", or something like that. See http://spatialreference.org/ for more info. If you aren't sure try this site: http://projfinder.com/ You can put the map on where the data is supposed to be, enter a coordinate pair, and it'll tell you what the best fit projection is for those coordinates.
## 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')
这篇关于从栅格中提取值时返回的NA的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
