CODEWITHSUNDEEP
rr-rasterr-sprgdal
Nigel Stackhouse
Nigel Stackhouse

Reputation: 492

Re-projecting a SpatialPointsDataFrame does not change the extent

I'm trying to overlay data related to bat locations (SpatialPointsDataFrame) on to the state of Colorado (SpatialPolygonsDataFrame). The CRS of the two objects are different:

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

When I reproject the bat data, the CRS does indeed change, but the extent is unaffected.

Before:

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

After:

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      : 139996.3, 748812, 4094998, 4535103  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#variables   : 17

I can't plot the points on the polygon, as the extent of my polygon object is dissimilar to that of my points object:

colorado
#class       : SpatialPolygonsDataFrame 
#features    : 1 
#extent      : -109.0608, -102.042, 36.99223, 41.00561  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#variables   : 13

Why is the extent of the SpatialPointsDataFrame not changing when re-projected and transformed?

Reproducible example(I tested this myself, and sent it to a friend; it IS reproducible with what is provided):

#Libraries
x = c('raster','sp','rgdal')

# If you don't know if you've installed the packages for some of these 
# libraries, run this:
# install.packages(x)

lapply(x, library, character.only=TRUE)
rm(x)

# Read in and format data -------------------------------------------------
BatRange_data = new("SpatialPoints", coords = structure(c(179935.719907205, 179935.938979813, 179935.938979813, 176598.335967664, 176598.335967664, 4499963.43180688, 4499963.30060606, 4499963.30060606, 4489332.4211975, 4489332.4211975), .Dim = c(5L, 2L), .Dimnames = list(NULL, c("coords.x1", "coords.x2"))), bbox = structure(c(176598.335967664, 4489332.4211975, 179935.938979813, 4499963.43180688), .Dim = c(2L, 2L), .Dimnames = list(c("coords.x1", "coords.x2"), c("min", "max"))) , proj4string = new("CRS"    , projargs = NA_character_))

#Use state bounds from gadm website:
us = getData("GADM", country="USA", level=1)

#Extract states (need to uppercase everything)
co = "Colorado" 

colorado = us[match(toupper(co),toupper(us$NAME_1)),]

#Re-project bat data to 'longlat'
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)

plot(BatRange_trnsfrmd)
plot(colorado,add=TRUE)

Upvotes: 2

Views: 6084

Answers (1)

lbusett
lbusett

Reputation: 5932

In your code:

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

you are overwriting the original projection of BatRange_data to lat/long before trying to reproject it. spTransform therefore is doing nothing, since it's trying to reproject from epsg:4326 to epsg:4326. This is why your extent is not changing.

You should therefore probably just remove this line:

crs(BatRange_data) = geo_proj

and all should work.

HTH.

Upvotes: 5

Related Questions