Projection problems with leaflet

[rlucas]

Hello,

I'm a R rookie trying to reproject a shape in order to display it with leaflet in a shiny application. My initial projection is

[+proj=utm +zone=22 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]

I tried several configs that I found on the web, but sadly I didn't found any infos on how to use the CRS() function

testproj1 <- spTransform(test,CRS('+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0'))

testproj2 <- spTransform(test, CRS("+proj=longlat +ellps=GRS80"))

testproj3 <- spTransform(test, CRS("+proj=longlat +datum=WGS84 +no_defs"))

testproj4 <- spTransform(test, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

Any of those gave a

Is projected: FALSE

when I check with a summary(spdf)

Anyone know how to define the CRS argument or wich one should I use?

[paleolimbot]

The reason that you are getting Is projected: FALSE is because all of the projections you are transforming to are longitude/latitude. Because longitude and latitude are angles (they are measured in degrees), they are not considered a "projected" coordinate system (even though your transform worked). This is what you want for use with leaflet! Except you want WGS84 latitude/longitude (your testproj3, which you can define as CRS("+init=epsg:4326") if you are using the sp package family. I would recommend using EPSG codes rather than learning PROJ syntax (what you have above).

library(sp)
shape_file <- system.file("shape/nc.shp", package = "sf")
spdf_original <- rgdal::readOGR(shape_file)
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/Library/Frameworks/R.framework/Versions/3.5/Resources/library/sf/shape/nc.shp", layer: "nc"
#> with 100 features
#> It has 14 fields
spdf_lon_lat <- spTransform(spdf_original, CRS("+init=epsg:4326"))

^{Created on 2019-04-04 by the reprex package (v0.2.1)}

I would also recommend using the sf family of packages rather than the sp, as it is much better documented (sf is also supported by leaflet). In the sf package family, you can use an EPSG code instead of a PROJ string to define a projection (in this case, 4326).

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
shape_file <- system.file("shape/nc.shp", package = "sf")
spdf_original <- read_sf(shape_file)
spdf_lon_lat <- st_transform(spdf_original, 4326)

^{Created on 2019-04-04 by the reprex package (v0.2.1)}

[rlucas]

Thanks for the tips. Its a good thing to know that long/lat are considered FALSE projection.
However it seem that it wasn't the reason of leaflet not displaying the data.

I use to connect to a PostGis database to get my data, using the custom function dbreadspatial from this website : https://neocarto.hypotheses.org/1186
This function directly connect to the database and import gis data into a spatialdataframe, according to his type (point, line, polygon...) and CRS.
So when I summarize my spdf, I get :

summary(test)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x 99404 384674
y 233785 493775
Is projected: TRUE
proj4string :
[+proj=utm +zone=22 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]

Then I apply the reprojection :

testproj <- spTransform(test, CRS("+proj=longlat +init=epsg:4326 +ellps=WGS84 +datum=WGS84 +no_defs"))

Wich looks good for use

summary(testproj)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -54.600647 -52.038341
y 2.111988 4.463896
Is projected: FALSE
proj4string :
[+proj=longlat +init=epsg:4326 +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0]

But when I try to display it on leaflet, I get an error message saying that I'm not using the long lat projection, even tought I do.

m <- leaflet() %>%

• addTiles() %>%
  

• setView(lng=-52.3333300, lat=4.9333300 , zoom=5) %>%
  

• addPolygons(map=m, data=sample_test, weight=2, col="black", opacity=0.5)
  

Error in derivePolygons(data, lng, lat, missing(lng), missing(lat), "addPolygons") :
addPolygons must be called with both lng and lat, or with neither.

So what's wrong with my method? I tried to connect to database with readOGR and applied the same transform and it didn't work either, with same error message.
I tried to reproject with st_transform but It didn't work :

sample_testproj <- st_transform(sample_test, 4326)
Error in UseMethod("st_transform") : no applicable method for 'st_transform' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialPolygonsNULL', 'SpatialVector')"

[paleolimbot]

In your first bit of code, I believe the problem is that you have map = m in your last line. Because you're using the pipe, the map argument is already supplied to addPolygons(), so you don't need to (actually you can't) supply the map argument yourself.

In your second bit of code, you're mixing sp/rgdal and sf. These are two different stacks of tools. To use sf, use read_sf() rather than readOGR(). As far as I remember, the arguments are the same for both functions, but read_sf() returns an object that works with the rest of the sf package.

[rlucas]

Thanks a lot, you were right on both sides. So the leaflet displaying problem was the map = m. And the transformation work with both sf or sp/rgdal packages if I do not mix it. To resume :

PostGIS database connexion

drv <- dbDriver("PostgreSQL")
dsn <- "PG:dbname='databasename' user = 'admin' password = 'pwd' host = '000.000.00.000' port = '5432'"

Importing/reprojecting data

sp/rgdal method

test <- readOGR(dsn=dsn, layer="layername")
testproj <- spTransform(test, CRS("+proj=longlat +init=epsg:4326 +ellps=WGS84 +datum=WGS84 +no_defs"))

sf method :

test <- read_sf(dsn=dsn, layer="layername")
testproj <- st_transform(test, 4326)

displaying on leaflet

m <- leaflet() %>%
addTiles() %>%
setView(lng=-52.3333300, lat=4.9333300 , zoom=5) %>%
addPolygons(data=sample_test, weight=2, col="black", opacity=0.5)
m

Thanks again!

[rlucas]

Last comment on the topic, I just found how to convert sp dataframe to sf dataframe on this excellent website wich detail all struggles that you can encounter working with tiddyverse and spatialdataframes.
https://geocompr.github.io/geocompkg/articles/tidyverse-pitfalls.html#pitfall-tidyverse-and-sp-dont-play

> library(spData)
> world_sp = as(world, "Spatial")
> world_sf = st_as_sf(world_sp)