0

I'm trying to plot point data (coördinates) on a multipolygoon map. But I can't get them on top of each other and in the right shape.

Underlying projection data of the Netherlands multipolygon data:

Simple feature collection with 1 feature and 3 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 619232.6
epsg (SRID):    28992
proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.87740,4.0725 +units=m +no_defs

Downloadlink to a .zip with 4 Geopackage files from the Netherlands

Underlying point projection data of coordinates:

Simple feature collection with 1466 features and 4 fields
    geometry type:  POINT
    dimension:      XY
    bbox:           xmin: 3.389737 ymin: 50.77106 xmax: 7.187327 ymax: 53.44402
    epsg (SRID):    28992
    proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.87740,4.0725 +units=m +no_defs

some rows of data:

points <- structure(list(brouwerijnaam = c("Almelosche", "HertogJan", "HertogJan", 
"Bavaria", "Bavaria", "Brand"), begin_datum_jaar = c(1987, 1981, 
1981, 1719, 1719, 1340), eind_datum_jaar = c(1990, 0, 0, 0, 0, 
0), plaats = c("Almelo", "Arcen", "Arcen", "Lieshout", "Lieshout", 
"Wijlre"), geometry = structure(list(structure(c(52.35192, 6.660562
), class = c("XY", "POINT", "sfg")), structure(c(51.473427, 6.181219
), class = c("XY", "POINT", "sfg")), structure(c(51.473427, 6.181219
), class = c("XY", "POINT", "sfg")), structure(c(51.518607, 5.597416
), class = c("XY", "POINT", "sfg")), structure(c(51.518607, 5.597416
), class = c("XY", "POINT", "sfg")), structure(c(50.832695, 5.895593
), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 50.832695, ymin = 5.597416, 
xmax = 52.35192, ymax = 6.660562), class = "bbox"), crs = structure(list(
    epsg = NA_integer_, proj4string = NA_character_), class = "crs"), n_empty = 0L)), row.names = c(NA, 
-6L), sf_column = "geometry", agr = structure(c(brouwerijnaam = NA_integer_, 
begin_datum_jaar = NA_integer_, eind_datum_jaar = NA_integer_, 
plaats = NA_integer_), .Label = c("constant", "aggregate", "identity"
), class = "factor"), class = c("sf", "tbl_df", "tbl", "data.frame"
))

Code I used to shape the x and y coördinates in a sf object with the same crs of multipolygon data.

library(sf)
library(tidyverse)

points <- df %>%
  filter(!is.na(lon)) %>%
  st_as_sf(coords = c("lat", "lon")) 

points <- st_set_crs(brouwerijen_sf, 28992)

ggplot() +
  geom_sf(data = multipolygon) +
  geom_sf(data = points) +
  coord_sf(datum = NA)

multipolygon Point data enter image description here

0

Without your data to be sure, here's my guess.

You need lon, then lat.

points <- df %>%
 filter(!is.na(lon)) %>%
 st_as_sf(coords = c("lon", "lat")) 

This is wrong:

points <- st_set_crs(brouwerijen_sf, 28992)

that assigned an Oblique Stereographic Alternative (+proj=sterea) projection to your multipolygon layer, and gave it the name "points".

You want to assign your points "longlat", 4326 is the EPSG shortcut.

points <- st_set_crs(points, 4326)

Now it should work. You can st_transform(points, st_crs(brouwerijen_sf)) in order to transform the points to the projection in use, but geom_sf will do that for you anyhow.

HTH

  • Sorry, made a mistake in my question brouwerijen_sf should be points. – Tdebeus Mar 25 at 11:56
  • setting 28992 on longlat is also wrong, can you update the Q? – mdsumner Mar 25 at 22:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.