rgisr-leaflet

Converting non-lat/lon coordinates into lat/lon coordinates in R

I would like to convert these non-lat/lon coords into lat/lon so I can map the on an interactive Leaflet map. I am not sure how to interpret the proj4string() output but I guess that is the key.

> a = readOGR("R03_11_WGS84","R03_11_WGS84")
OGR data source with driver: ESRI Shapefile 
Source: "R03_11_WGS84", layer: "R03_11_WGS84"
with 53173 features
It has 24 fields
Integer64 fields read as strings:  SEZ2011 
> a1 = a@polygons[[1]]@Polygons[[1]]@coords
> head(a1)
         [,1]    [,2]
[1,] 486771.9 4999855
[2,] 486807.4 4999833
[3,] 486825.8 4999826
[4,] 486843.7 4999821
[5,] 486856.2 4999818
[6,] 486884.5 4999810
> proj4string(a)
[1] "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

The blue line on the top is the polygon I added from the SPDF:

enter image description here

From what I read there could be a way to change leaflet's CRS.

Solution

  • Use spTransform. The required proj4string for using lat/lon in Leaflet is "+proj=longlat +datum=WGS84"

    library(rgdal)
spTransform(a, CRS("+proj=longlat +datum=WGS84"))

    If that doesn't work, the official definition of EPSG:3857 is

    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"

    Reprex

    url <- "http://www.istat.it/storage/cartografia/basi_territoriali/WGS_84_UTM/2011/R03_11_WGS84.zip"
destfile <- basename(url)
destfolder <- tools::file_path_sans_ext(destfile)
download.file(url = url, destfile = destfile)
unzip(destfile)

library(rgdal)
#> Loading required package: sp
#> rgdal: version: 1.2-15, (SVN revision 691)
#>  Geospatial Data Abstraction Library extensions to R successfully loaded
#>  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
#>  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
#>  GDAL binary built with GEOS: FALSE 
#>  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
#>  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
#>  Linking to sp version: 1.2-5
a <-  readOGR(dsn = destfolder, layer = destfolder)
#> OGR data source with driver: ESRI Shapefile 
#> Source: "R03_11_WGS84", layer: "R03_11_WGS84"
#> with 53173 features
#> It has 24 fields
#> Integer64 fields read as strings:  SEZ2011
a_coords <- a@polygons[[1]]@Polygons[[1]]@coords
head(a_coords)
#>          [,1]    [,2]
#> [1,] 486771.9 4999855
#> [2,] 486807.4 4999833
#> [3,] 486825.8 4999826
#> [4,] 486843.7 4999821
#> [5,] 486856.2 4999818
#> [6,] 486884.5 4999810

b <- spTransform(a, CRS("+proj=longlat +datum=WGS84"))
b_coords <- b@polygons[[1]]@Polygons[[1]]@coords
head(b_coords)
#>          [,1]     [,2]
#> [1,] 8.831718 45.15205
#> [2,] 8.832170 45.15185
#> [3,] 8.832404 45.15179
#> [4,] 8.832632 45.15174
#> [5,] 8.832791 45.15171
#> [6,] 8.833151 45.15164