Search code examples
rleafletr-sprgdalr-leaflet

Misalignment between POSGAR94 polygons and WGS84 leaflet map


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:

actual map

All the polygons are offset by a block, its mostly visible on the city perimeter. Here's how that polygon should look like:

expected map

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"]]

Solution

  • 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:

    BSAS

    But going a bit south (e.g. around Mar del Plata coastline) shows an obvious misalignment:

    Mar del Plata

    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:

    ideign + indec

    In fact, superimposing OSM, IDE-IGN and INDEC data means three different data sources which don't match:

    ideign + indec + osm

    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.