0

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"]]
  • I recently had similar issues with my basemap and coordinates not matching up and this was the trick: https://github.com/dkahle/ggmap/issues/160 I realize this is ggmap, not leaftet, but I hope it helps. – TTS Jul 09 '20 at 14:34
  • I'll play around with that idea, but it seems that just google had that problem while using ggmap:get_map() instead of dismo:gmap(). Maybe its the same with OpenStreetMap.. – Catriel Lopez Jul 09 '20 at 15:02
  • The [dataset you linked to](https://www.indec.gob.ar/ftp/cuadros/territorio/codgeo/Codgeo_Pais_x_prov_datos.zip) is a shapefile in EPSG:4326 with data for provinces, but the problem (as indicated by your screenshots) is about postal codes (or neighbourhoods? districts?) in EPSG:22183. Please edit your question and point to the problematic dataset. – IvanSanchez Jul 09 '20 at 15:03
  • Edited! [This is the correct dataset](https://www.indec.gob.ar/ftp/cuadros/territorio/codgeo/Codgeo_Buenos_Aires_con_datos.zip) – Catriel Lopez Jul 09 '20 at 15:21

1 Answers1

0

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.

IvanSanchez
  • 18,272
  • 3
  • 30
  • 45