I need to draw a bunch of polygons from this dataset on a leaflet map on R:
The coordinates are in POSGAR94, but I need them in WGS84 to plot on leaflet (over a OpenStreetMap layer) and to compare them with other data (already on WGS84):
library(rgdal)
library(magrittr)
library(leaflet)
complete_data <- readOGR("data_folder",
GDAL1_integer64_policy = TRUE)
complete_data <- spTransform(bsas,
CRS("+proj=longlat +datum=WGS84 +no_defs"))
I filter the data to keep only a section of it:
int_data <- complete_data[grep("^0604219|^0604201|060421102|060421103", complete_data@data$link), ]
And I plot:
leaflet(int_data, options = leafletOptions(minZoom = 12, maxZoom = 18)) %>%
setMaxBounds(lat1 = -37.1815, lng1 = -58.5581, lat2 = -37.1197, lng2 = -58.4297) %>%
addTiles()%>%
addPolygons(color = "#3498db", weight = 2, smoothFactor = 0.5,
opacity = 0.5, fillOpacity = 0.1,
highlightOptions = highlightOptions(color = "black", weight = 3,
bringToFront = TRUE))
The current result looks like:
All the polygons are offset by a block, its mostly visible on the city perimeter. Here's how that polygon should look like:
My questions are:
Am I making a mistake with the proyection? Or does spTransform introduce an error in the coordinates?
or
Is my code ok, but the data is wrong?
EDIT: This is the output of st_crs before and after the conversion:
BEFORE
Coordinate Reference System:
User input: +proj=tmerc +lat_0=-90 +lon_0=-66 +k=1 +x_0=3500000 +y_0=0 +ellps=WGS84 +units=m +no_defs
wkt:
PROJCS["unnamed",
GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",-90],
PARAMETER["central_meridian",-66],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",3500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
AFTER
Coordinate Reference System:
User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
wkt:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]
This looks like an issue with the dataset. I'm looking at it in QGIS, along with some OSM basemap, and around Buenos Aires everything seems to fit the road network nicely:
But going a bit south (e.g. around Mar del Plata coastline) shows an obvious misalignment:
Is my code ok, but the data is wrong?
Since the same problem can be seen using a completely different method, it's safe to say that your code is OK, and working as expected.
We can say that the dataset mismatches when reprojected over a OSM basemap. However, we cannot say that the data is wrong. If we load your dataset along with some reference data from ide.ign.gob.ar (department limits), the data also doesn't match:
In fact, superimposing OSM, IDE-IGN and INDEC data means three different data sources which don't match:
This is normal. There is no easy definition of "right" when it comes to alignment of GIS datasets, and there are a lot of factors in play: collection criteria, orthorectification parameters, datum shifts, projection warping and even continental drift, among others.