Search code examples
rgeolocationgeospatialr-sf

Transform point into UTM 32N

i am coding with R and RStudio. I have a set of polygons in one in my global environment called buek250_shapefile. I am able to plot them with leaflet. For the next step I wrote a function which gets lat and lon values, and should give back information about the shape which haqs the shortest distance to the next polygon shape. My problem: I always get back the same polygon shape for all lat lon values I am trying. And the computed distances are between 5000 and 7000km. But all my shapes and lat lon values are in Germany, and Germany is not so big. My data: https://download.bgr.de/bgr/boden/BUEK250/shp/buek250_mg_utm_v55.zip Can anyone here help me please? Thank you very much

lowestdistancebuek <- function(inputlat, inputlon) {
  # Modify the clicked point. Make the input data into CRS format
  # Order is important: first lat, then lon
  point <- st_point(c(inputlat, inputlon)) %>%
    st_sfc(crs = 25832)
  
  # Is there a shape containing the clickedpoint?
  nearestshape <- st_intersects(point, buek250_shapefile_all, sparse = FALSE)
  # nearestsape is a logical matrix. For every polygon it contains if the point is in or not
  # Which polygon contains the point? Extract it
  nearestshape <- which(nearestshape, arr.ind = TRUE)
  # Extract the column
  nearestshape <- nearestshape[,2]

  # If there is a shape containing the point, nearestshape has length 1
  # In this case, nearestshape is set to the sf
  # Otherwise, the nearestshape is determinded and set to nearestshape
  if(length(nearestshape) == 1) {
    nearestshape <- buek250_shapefile_all[nearestshape, ]
  } else if (length(nearestshape) == 0){ # no polygon contains the point
    # Compute the distances between each polygon and the point. Result is saved in distances
    distances <- st_distance(point, buek250_shapefile_all,
                             by_element = FALSE, # TRUE returns a vector
                             tolerance = TRUE # speeds up computing
                             )
    # Which distance ist the smallest one?
    nearestshape <- which.min(distances)
    # The shortest number was saved in nearestshape, now there is saved the polygone information
    nearestshape <- buek250_shapefile_all[nearestshape, ]
  }
  return(nearestshape)
}

Right now I receive for every input the same output. I expect the code to return the shapefile with the shortest distance to the given point.

Solution

  • While geographic coordinates are usually represented as latitude, longitude pairs, geometric operations generally expect x, y order, that also applies for sf::st_point(). If you need a mnemonic for lat/lon, you could try "lat" -> "ladder" -> "y" .

    To transform from one CRS to another, use st_transform(); but for this to work, input must have CRS assigned and this is the role of crs argument in st_sfc(x, crs = "WGS84") .It defines what you currently have, not what you want to end up with.

    I'm using shapes from giscoR for a reproducible example here.

    library(giscoR)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

# example shapes from GISCO / giscoR
ger_25832 <- 
  giscoR::gisco_get_nuts(country = "Germany", nuts_level = "1") |>
  st_transform(25832)

plot(ger_25832[,"NUTS_NAME"], axes = TRUE)

    First, let's try to build a point from those geographic coordinates so it ends up somewhere in the area of interest:

    # input point coordinates, WGS84
inputlat_wgs84 <- 52.52   # mnemonic: latitude like "ladder", y-axis
inputlon_wgs84 <- 10.405

# create WGS84 point, transform, check distances:
ref_pt <- 
  # x,y!
  st_point(c(inputlon_wgs84, inputlat_wgs84)) |> 
  st_sfc(crs = "WGS84") |>
  st_transform(crs = st_crs(ger_25832))

st_distance(ger_25832, ref_pt) |> units::set_units(km)
#> Units: [km]
#>            [,1]
#>  [1,] 221.48755
#>  [2,] 184.04861
#>  [3,]  88.60661
#>  [4,] 114.36606
#>  [5,] 101.80472
#>  [6,] 115.41499
#>  [7,]  87.67142
#>  [8,]   0.00000
#>  [9,]  87.76893
#> [10,] 250.06811
#> [11,] 396.54742
#> [12,] 165.09594
#> [13,]  35.32782
#> [14,]  94.75578
#> [15,]  99.73750
#> [16,] 308.01952

    Seems that we got lat/lon order and transformation right this time, let's adjust that function accordingly, we can use st_nearest_feature() to subset nearest feature:

    # same in the function + subset with st_nearest_feature()
lowestdistancebuek <- function(inputlat_wgs84, inputlon_wgs84, shapes = ger_25832){
  ref_pt <- 
    st_point(c(inputlon_wgs84, inputlat_wgs84)) |> 
    st_sfc(crs = "WGS84") |>
    st_transform(crs = st_crs(shapes))
  shapes[st_nearest_feature(ref_pt, shapes),]
}

nearest_shape <- lowestdistancebuek(inputlat_wgs84, inputlon_wgs84)
nearest_shape
#> Simple feature collection with 1 feature and 10 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 343686.6 ymin: 5682911 xmax: 674176.7 ymax: 5971541
#> Projected CRS: ETRS89 / UTM zone 32N
#>    NUTS_ID LEVL_CODE URBN_TYPE CNTR_CODE     NAME_LATN     NUTS_NAME MOUNT_TYPE
#> 58     DE9         1         0        DE NIEDERSACHSEN NIEDERSACHSEN          0
#>    COAST_TYPE FID geo                       geometry
#> 58          0 DE9 DE9 MULTIPOLYGON (((486772.8 59...

    Plot our test point & resulting shape:

    # plot 
plot(st_union(ger_25832))
plot(nearest_shape, add = TRUE)
#> Warning in plot.sf(nearest_shape, add = TRUE): ignoring all but the first
#> attribute
plot(ref_pt, col = "red", pch = 16, add = TRUE)

    Created on 2024-06-19 with reprex v2.1.0