2017-10-12 6 views
-1

J'essaie de superposer des données liées à l'emplacement des chauves-souris (SpatialPointsDataFrame) sur l'état du Colorado (SpatialPolygonsDataFrame). Le CRS des deux objets sont différents:Re-projeter un SpatialPointsDataFrame ne change pas l'étendue

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 

Quand je reprojeter les données de chauve-souris, le CRS ne change en effet, mais la mesure n'est pas affecté.

Avant:

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 

Après:

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 

Je ne peux pas tracer les points sur le polygone, comme l'étendue de mon objet polygone est très différent de celui de mon objet de points:

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 

Pourquoi l'étendue de la SpatialPointsDataFrame ne change-t-elle pas lorsqu'elle est re-projetée et transformée?

exemple Reproductibles (J'ai testé moi-même, et envoyé à un ami, il est reproductible avec ce qui est prévu):

#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) 
+0

L'exemple n'est pas reproductible pour nous, car il dépend de "Armstrong_Edited"; et chargement inutilement compliqué d'environ 30 bibliothèques, alors que trois suffiraient. (ou deux: 'library (raster); library (rgdal)') – RobertH

+0

"Armstrong_Edited" est utilisé pour faire 'BatRange_data', pour lequel j'inclus la sortie' dput' à la fin. Donc, ceci EST reproductible pour vous. Si vous ne souhaitez pas charger toutes les bibliothèques, supprimez celles que vous ne voulez pas. Cela prend quelques secondes à charger, donc ce n'est pas vraiment "compliqué". –

Répondre

2

Dans votre 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 

vous êtes écrasant la projection originale de BatRange_data à lat/long avant en essayant de le reprojeter. spTransform ne fait donc rien, puisqu'il essaie de reprojeter epsg: 4326 à epsg: 4326. C'est pourquoi votre étendue ne change pas.

Vous devriez donc probablement juste supprimer cette ligne:

crs(BatRange_data) = geo_proj

et tout devrait fonctionner.

HTH.

+0

Merci, cela a fait le tour exactement! –