'Extracting Raster Variables but only getting NA's
I am trying to extract variable values from 16 .img files and bind it to a pre-existing structure, and while the code runs, it only spits back NA values. I've tried converting the crs but it hasn't done anything.
library(raster)
library(sp)
library(rgdal)
library(sf)
library(exactextractr)
setwd("~/preds")
pres.abs <- st_read("~spp106mort.PPsA.SF.shp") # load in GIS data
etpt_5.r <- raster("etpt_5.img", band=1)
etpt_6.r <- raster("etpt_6.img", band=1)
etpt_sprin.r <- raster("etpt_sprin.img", band=1)
exp1nrm.r <- raster("exp1nrm.img", band=1)
exp3nrm.r <- raster("exp3nrm.img", band=1)
exp5nrm.r <- raster("exp5nrm.img", band=1)
mind_yr_av.r <- raster("mind_yr_av.img", band=1)
prad_sw_di.r <- raster("prad_sw_di.img", band=1)
prec_w_hal.r <- raster("prec_w_hal.img", band=1)
prec_winte.r <- raster("prec_winte.img", band=1)
rough_1k.r <- raster("rough_1k.img", band=1)
tave_s_hal.r <- raster("tave_s_hal.img", band=1)
tave_sprin.r <- raster("tave_sprin.img", band=1)
tmax_s_hal.r <- raster("tmax_s_hal.img", band=1)
tmax_summe.r <- raster("tmax_summe.img", band=1)
topos.r <- raster("topos.img", band=1)
allvars <- stack(etpt_5.r, etpt_6.r, etpt_sprin.r, exp1nrm.r, exp3nrm.r, exp5nrm.r,
mind_yr_av.r, prad_sw_di.r, prec_w_hal.r, prec_winte.r, rough_1k.r,
tave_s_hal.r, tave_sprin.r, tmax_s_hal.r, tmax_summe.r, topos.r)
## converting pres.abs to new crs
pres.abs.old <- pres.abs
newcrs <- crs(allvars)
pres.abs.new <- st_transform(pres.abs, newcrs)
crs(allvars)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["WGS 84 (with axis order normalized for visualization)",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic longitude (Lon)",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["geodetic latitude (Lat)",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
projection(allvars)
[1] "+proj=longlat +datum=WGS84 +no_defs"
crs(pres.abs.new)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
WKT2 2019 representation:
GEOGCRS["WGS 84 (with axis order normalized for visualization)",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic longitude (Lon)",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["geodetic latitude (Lat)",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
projection(pres.abs.new)
[1] "+proj=longlat +datum=WGS84 +no_defs"
So I got the crs and projection the same, but when it comes time to extract, here's what happens
head(pres.abs.new)
Simple feature collection with 6 features and 12 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -96.00105 ymin: 23.00031 xmax: -96.00105 ymax: 23.00031
Geodetic CRS: WGS 84 (with axis order normalized for visualization)
UNIQUEID wgs_xF wgs_yF MORT106 exp5nrm topos rough_1k prad_sw_di tmax_s_hal etpt_6 mind_yr_av prec_w_hal
1 Z0401030100988520 -109.6594 33.25329 1 -23 -1.307692 25.61479 -5511 262.3333 101 -1676.667 34.83333
2 Z0401030301183589 -109.4448 33.20588 1 -48 -5.153846 51.95938 -5622 278.6667 101 -1708.083 29.83333
3 Z0401030401180927 -109.2798 33.42938 1 21 19.923077 39.72136 -5198 258.0000 100 -1680.667 33.50000
4 Z0401030501183471 -109.3880 33.38541 1 0 18.230770 27.97424 -5506 250.8333 100 -1585.667 35.83333
5 Z0401030601186305 -109.3887 33.34027 1 9 24.153847 34.67730 -5456 262.8333 101 -1679.917 33.00000
6 Z0401030701187126 -109.2232 33.56368 1 110 -33.538460 61.15303 -5409 248.5000 99 -1495.417 34.16667
geometry
1 POINT (-96.00105 23.00031)
2 POINT (-96.00105 23.00031)
3 POINT (-96.00105 23.00031)
4 POINT (-96.00105 23.00031)
5 POINT (-96.00105 23.00031)
6 POINT (-96.00105 23.00031)
extracted <- extract(allvars, pres.abs.new[, c("wgs_xF", "wgs_yF")])
etpt_5 etpt_6 etpt_sprin exp1nrm exp3nrm exp5nrm mind_yr_av prad_sw_di prec_w_hal prec_winte rough_1k tave_s_hal tave_sprin tmax_s_hal
[1,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[2,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[3,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[4,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[5,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[6,] NA NA NA NA NA NA NA NA NA NA NA NA NA NA
tmax_summe topos
[1,] NA NA
[2,] NA NA
[3,] NA NA
[4,] NA NA
[5,] NA NA
[6,] NA NA
This is weird because each of the rasters I imported have distinct minimum and maximum values, so I know that they're not actually NA
prec_winte.r
class : RasterLayer
dimensions : 1518, 1478, 2243604 (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333 (x, y)
extent : -114.6167, -102.3, 29.49167, 42.14167 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : prec_winte.img
names : prec_winte
values : 3.333333, 140.3333 (min, max)
I think the issue is dependent on where in the raster or dataframe I'm pulling the values from, but I can't figure out where!
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|