Search code examples
rgisr-sfgeosphere

Generating GPS Coordinates / creating a raster of evenly distributed points 1km apart


i would like to cover an area defined by a bbox lat long coordinates with a raster of gps points 1 km apart. Currently i generate 2000 points for a bboxbbox=8.9771580802,47.2703623267,13.8350427083,50.5644529365 the following way:

as.data.frame(cbind(runif(2000,8.9771580802 ,13.8350427083),runif(2000,47.2703623267,50.5644529365)))

Since runif is a normal distribution, i think i just have to increase the amount of points to cover the whole area the way i need it. Is there a more clever way to do it? How many points would i need?

UPDATE

I thought i maybe can use the package sp to do the job but im still not realy familiar the with settings:

  longitudes <- c(8.9771580802, 13.8350427083)
  latitudes <- c(47.2703623267, 50.5644529365)
  bounding_box <- matrix(c(longitudes, latitudes), nrow = 2, byrow = TRUE, dimnames = list(NULL, c("min", "max")))
  projection <- "+proj=longlat" 
  sp_box<-Spatial(bbox = bounding_box, proj4string = CRS(projection))
  p_sample<-spsample(sp_box, 10000, type="regular")

If i understand correctly this will give me a number of points evenly distributed within my coordinates. spsample has an option for cell size but i dont grasp it yet. BR Andreas


Solution

  • I'm not too much into spatial data and analysis but maybe it helps as a first step (I took different coordinates to get a reproducible example, which fits the dimensions of Germany where I have some feelings for the dimensions). I'm sure that there is a more elegant way but it should give you what you need. geosphere::destPoint() is ued to compute the points for a given distance and direction, geosphere::distGeo() computes the north-south/west-east distance of the given box to compute how many points we need for each direction. expand_grid() is then used to compute every combination for the computed border points.

    Be also aware that I changed the distance between the points to 10,000 meters or 10 km to get fewer points and a nicer plot. You would have to change the numbers accordingly

    nw <- c(5.8 55) 
    se <- c(15.1, 47)
    
    lon1 <- nw[1]
    lat1 <- nw[2]
    lon2 <- se[1]
    lat2 <- se[2]
    
    
    #(1) compute the border points in y direction, going south from the nw-point 
    # while keeping lon constant
    
    lat <- geosphere::destPoint(nw, 180, 1:floor(geosphere::distGeo(c(lon1,lat1), 
                                                                     c(lon1,lat2))/10000)*10000) 
    lat <- as_tibble(lat) 
    
    #(2) compute the border point in x direction (analog to above)
    lon <- geosphere::destPoint(nw, 90, 1:floor(geosphere::distGeo(c(lon1,lat1), 
                                                                    c(lon2,lat1))/10000)*10000)
    lon <- as_tibble(lon)
    
    # use expand_grid() to compute all combinations 
    grid <- tidyr::expand_grid(lat$lat, lon$lon) 
    names(grid) <- c("lat", "lon") #nicer names
    
    ### for visualizing what we've done, map germany with a grid overlay
    germany  <- rnaturalearth::ne_countries(type =  "countries", 
                country = "germany", returnclass = "sf")
    
    ggplot2::ggplot(data = germany)+
      ggplot2::geom_sf()+
      ggplot2::geom_point(data = grid, mapping = aes(x = lon, y = lat), size = 0.01)
    

    Map of Germany overlaid with equidistant grid