CODEWITHSUNDEEP
rcsvgisshapefile
PVic
PVic

Reputation: 449

R Map CSV & SHP

I'm trying to plot latitude & longitude coordinates on top of a shapefile. I am told they have the same coordinate system (NAD27). I have tried numerous ways after researching it, but nothing works perfectly. They either map one or the other, but not both. Below is an example of reading them in, but I don't even know how to begin starting to plot since my csv is from the sp package my the shp is from the maptools package.

pts <- read.csv("locations.csv")
pts <- pts[,4:5]
names(pts) <- c("Latitude", "Longitude")
pts <- pts[complete.cases(pts),]
pts <- SpatialPoints(pts)
border <- maptools::readShapePoly("landman_one shape.shp")

Any help would be greatly appreciated.

Below are the different structures.

> str(pts)
Formal class 'SpatialPoints' [package "sp"] with 3 slots
  ..@ coords     : num [1:372, 1:2] 30.9 30.9 31 31 31 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:372] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "Latitude" "Longitude"
  ..@ bbox       : num [1:2, 1:2] 29.5 -105.2 36 -76.1
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "Latitude" "Longitude"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA

> str(border)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1 obs. of  2 variables:
  .. ..$ Id       : int 0
  .. ..$ PERIM_GEO: num 998432
  .. ..- attr(*, "data_types")= chr [1:2] "N" "F"
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 922511 545445
  .. .. .. .. .. .. ..@ area   : num 3.45e+10
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:252, 1:2] 954182 954073 954006 953914 953828 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 922511 545445
  .. .. .. ..@ ID       : chr "0"
  .. .. .. ..@ area     : num 3.45e+10
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 819851 386743 1042230 669865
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA

Also, I'm not too familiar with projections, so any guidance with that would be appreciated. Also, if there's a way to do this with the same package I'd be open to that.

Upvotes: 0

Views: 177

Answers (1)

Phil
Phil

Reputation: 4444

The problem is that neither of your spatial data sets have a projection set and are using different coordinate systems. You can see this in the @proj4string line, which are both currently NA. A quick way to ascertain that the coordinate systems are different is with the @bbox slot ('bounding box'). In pts these are longitude/latitude as they are within -180/+180 and -90/+90 respectively. border, on the other hand, has values that are in entirely different units (819851 386743 1042230 669865). Neither of these look to be in the NAD27 coordinate system.

To solve this, first I would use rgdal::readOGR() or sf::read_sf() to read in your border shape file. These functions load in a projection if one is present, so this should tell you what projection this data is in.

Second, to set the projection for pts you can use proj4string():

proj4string(pts) = CRS("+init=epsg:4326")

Once both pts and border have a projection set, you can transform one to the other coordinate system with spTransform(). As I don't know what the CRS/projection for border is I'll convert this to WGS84 (the projection of pts), but just swap them to convert pts to the same projection as border once you know its CRS.

border = spTransform(border, CRS("+init=epsg:4326"))

Once in the same projection they should plot together. You can confirm their projections match with:

proj4string(border) == proj4string(pts)

Upvotes: 2

Related Questions