[GIS] Plot points by longitude and latitude on a map created by centroids

centroidsgeolocationlatitude longitudershapefile

I'm searching for a way to plot some points for which I have the longitude and latitude (using Google Maps API) on a map created with centroids. I succeeded to do it by using functions such as map_data() and functions from the package ggplot2.

Results (just to show it really works) :

However I don't like the french departments datas given by map_data(), it's quite ambiguous but the nice thing is that their coordinates are given by longitude and latitude. (french departments are similar to counties)

What I would like to do is to use a shapefile from a site such as GEOFLA (here) that provides me datas in projection LAMBERT-93 to plot my map. However the shapefile only contains the coordinates of deparments centroids so, since I have, in one hand, coordinates given by longitude and latitude and, in the other hand, coordinates given by centroids, so I can't plot the points on my map since the coordinates are not in the same "unit".

That's what I get with summary() on the shapefile :

Object of class SpatialPolygonsDataFrame
Coordinates:
        min     max
x   99217.1 1242417
y 6049646.3 7110480
Is projected: TRUE 
proj4string :
[+proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m
+no_defs]

Summary of map_data() (I show only the coordinates) :

      long               lat     
 Min.   :-5.14209   Min.   :41.37
 1st Qu.: 0.05213   1st Qu.:44.41
 Median : 2.45500   Median :46.79
 Mean   : 2.45851   Mean   :46.47
 3rd Qu.: 4.97283   3rd Qu.:48.44
 Max.   : 9.56267   Max.   :51.10

This shows that the coordinates used are different.

I tried many things such as using proj4string(), sptransform() for example to get the same units but I didn't get what I wanted.

So how should I do to get the same units for my two differents types of coordinates ? Is there a way to get the longitude and latitude from the centroids of the shapefile or should I use a function (if it existed) that can transform my points with longitude and latitude into centroids ?

This is my first time using R for maps and polygons, projections etc. also so don't be surprised if I said some real stupid things.

Exemple of points I want to plot :

  centerid latitude longitude                                       adresse_complete
1   0121H1 46.22368  5.210819               900 Route de Paris, 01440 Viriat, France
2   0121H2 46.22368  5.210819               900 Route de Paris, 01440 Viriat, France
3   0218H3 49.37008  3.335127 46 Avenue du Général de Gaulle, 02200 Soissons, France
4   0303H1 46.34292  2.608996       18 Avenue du 8 Mai 1945, 03100 Montluçon, France
5   0303H2 46.34292  2.608996       18 Avenue du 8 Mai 1945, 03100 Montluçon, France

dput() from above :

structure(list(centerid = c("0121H1", "0121H2", "0218H3", "0303H1", 
"0303H2"), latitude = c(46.2236804, 46.2236804, 49.3700842, 46.3429172, 
46.3429172), longitude = c(5.2108193, 5.2108193, 3.335127, 2.6089958, 
2.6089958), adresse_complete = c("900 Route de Paris, 01440 Viriat, France", 
"900 Route de Paris, 01440 Viriat, France", "46 Avenue du Général de Gaulle, 02200 Soissons, France", 
"18 Avenue du 8 Mai 1945, 03100 Montluçon, France", "18 Avenue du 8 Mai 1945, 03100 Montluçon, France"
)), .Names = c("centerid", "latitude", "longitude", "adresse_complete"
), row.names = c(NA, 5L), class = "data.frame")

Best Answer

You can either reproject your polygon shapefile OR your dataset. EPSG is a short code for spatial projections. EPSG:4326 is the code for non-projected data in degrees using wgs84: http://spatialreference.org/ref/epsg/4326/

Reproject your polygon into 4326

france_wgs84 <- spTransform(france_L93, CRS("+init=epsg:4326"))

Reproject your dataset into Lambert93

data <- structure(list(centerid = c("0121H1", "0121H2", "0218H3", "0303H1", 
"0303H2"), latitude = c(46.2236804, 46.2236804, 49.3700842, 46.3429172, 
46.3429172), longitude = c(5.2108193, 5.2108193, 3.335127, 2.6089958, 
2.6089958), adresse_complete = c("900 Route de Paris, 01440 Viriat, France", 
"900 Route de Paris, 01440 Viriat, France", "46 Avenue du Général de Gaulle, 02200 Soissons, France", 
"18 Avenue du 8 Mai 1945, 03100 Montluçon, France", "18 Avenue du 8 Mai 1945, 03100 Montluçon, France"
)), .Names = c("centerid", "latitude", "longitude", "adresse_complete"
), row.names = c(NA, 5L), class = "data.frame")


coordinates(data) <- ~longitude+latitude
proj4string(data) <- "+init=epsg:4326"
data_L93 <- spTransform(data, CRS("+proj=lcc +lat_1=44 +lat_2=49 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m
+no_defs"))

Get example data:

# Load libraries
library('rgdal')
library('raster')

# Get SpatialPolygonsDataFrame object example
polygon <- getData('GADM', country = 'URY', level = 0)

# Plot data
plot(polygon)

# Make sample points 
sample.spdf <- spsample(polygon, n = 200, type = "random")
sample.spdf <- spsample(sample.spdf, n = 200, type = "random") # get points outside polygon
sample.spdf <- SpatialPointsDataFrame(sample.spdf, data = data.frame("ID" = 1:200, "V2" = sample(x = 1:1000, size = 200))) # add random data

Plot data:

# Plot points and polygon
plot(polygon)
plot(sample.spdf, add = TRUE, pch = 19, cex = 1)
box()

# Interactive map with Mapview
library(mapview)
mapview(polygon) + mapview(sample.spdf)

Spatial Query: points in polygon with different CRS:

# Apply transformation to sample.spdf
sample.spdf.t <- spTransform(sample.spdf, CRSobj = "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Ask layers CRS
projection(sample.spdf.t) # [1] "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(polygon) # [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

pt.in.poly.sp <- sp::over(sample.spdf.t, polygon, fn = NULL) 
# Error: identicalCRS(x, y) is not TRUE
pt.in.poly.rgeos <- rgeos::gIntersection(sample.spdf.t, polygon) 
# spgeom1 and spgeom2 have different proj4 strings

# Plot points transformed + polygon
plot(polygon)
plot(sampleSP.t, add = TRUE, pch = 19, cex = 0.5)

Spatial Query: points in polygon with same CRS:

# Ask layers CRS
projection(sample.spdf) # [1] [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(polygon) # [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

pt.in.poly.sp <- sp::over(sample.spdf, polygon, fn = NULL)
pt.in.poly.rgeos <- rgeos::gIntersection(sample.spdf, polygon)

# Plot points in polygon
plot(polygon)
plot(sample.spdf, pch = 19, cex = 1, add = TRUE)
plot(pt.in.poly.rgeos, pch = 19, cex = 0.5, col = 'green', add = TRUE)
box()

# Interactive map with Mapview
mapview(polygon) + mapview(sample.spdf) + mapview(pt.in.poly.rgeos)

Shapefile in R – How to Calculate Area of Shapefile Polygons in Square Kilometers

An error in a factor of about ten in area corresponds to a scale error of about 3.16 in length (square root of ten), which I have seen happen when people use imperial foot units instead of metres (or yards for really old projection systems).

Oh look... https://epsg.io/2314

Remarks: Foot version of Trinidad 1903 / Trinidad Grid (code 30200) used by some US-based companies including Amoco Trinidad.

And your expanded EPSG string even has:

+to_meter=0.3047972654

(kicking myself for not noticing earlier...)

the projection system is in feet is not in square metres, so the area is in square feet. Scale it by feet^2 to meters^2. Your area is:

> 3236.118 * .3047972654^2
[1] 300.6398

Ta-da!

A bit more digging shows that the unit is "Clarke's Foot" which is a bit different to the International Foot, the US Surveyor's Foot, and the Indian Foot, which are the only feet supported as units by my current version of PROJ, hence the explicit scaling factor in the PROJ string.

Note that if you use sf spatial objects instead of sp, then once you've set the CRS it will compute areas and give the units of area, but make sure you are using the latest sf package since Edzer only fixed a bug that ignores +to_meters parameters today (after I reported it because I was trying sf with your projection).