CODEWITHSUNDEEP
rrandomgeospatialdata-manipulationshapefile
stats_noob
stats_noob

Reputation: 5897

Simulating Random Points in Different ZIP Codes

I have the following shapefile in R and created this map of eastern United States.

library(sf)  
library(leaflet)
library(leafgl)
library(colourvalues)
library(leaflet.extras)


nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>% 
  st_transform(st_crs(4326)) %>% 
  st_cast('POLYGON')

The shapefile looks something like this:

> nc
Simple feature collection with 108 features and 14 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.45662 ymax: 36.58973
Geodetic CRS:  WGS 84
First 10 features:
     AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79                           geom
1   0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1      10  1364     0      19 POLYGON ((-81.47258 36.2344...
2   0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0      10   542     3      12 POLYGON ((-81.23971 36.3654...
3   0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5     208  3616     6     260 POLYGON ((-80.45614 36.2426...
4   0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145 POLYGON ((-76.00863 36.3196...
4.1 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145 POLYGON ((-76.02682 36.5567...
4.2 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1     123   830     2     145 POLYGON ((-75.90164 36.5562...
5   0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9    1066  1606     3    1197 POLYGON ((-77.21736 36.2410...
6   0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7     954  1838     5    1237 POLYGON ((-76.74474 36.2339...
7   0.062     1.547  1834    1834      Camden 37029  37029       15   286     0     115   350     2     139 POLYGON ((-76.00863 36.3196...
8   0.091     1.284  1835    1835       Gates 37073  37073       37   420     0     254   594     2     371 POLYGON ((-76.56218 36.3406...

I would like to simulate a random longitude/latitude point that falls within the geographical confines of "Ashe" - but I am not sure how to do this.

I see that in this shapefile, there is a column called "geom" which appears to contain information on the geographical boundaries of each location (e.g. Ashe, Alleghany, Surry, etc.).

But is there something I can do to simulate random longitude/latitude points and determine where they are situated?

Conceptually, I thought of two options to accomplish this:

Option 1: As an example, suppose if I simulate the following points:

id = 1:100
latitude = rnorm(100,-81, 0.15)
longitude = rnorm(100,36.2, 0.15)
my_data = data.frame(id, latitude, longitude)

  id  latitude longitude
1  1 -81.15816  36.42389
2  2 -81.40090  36.23823
3  3 -80.97732  35.97633
4  4 -80.80150  36.20300
5  5 -81.26429  36.23899
6  6 -81.13721  36.31100

I would like to find out which areas these points are located in - for example:

  id  latitude longitude  location
1  1 -81.15816  36.42389      Ashe
2  2 -81.40090  36.23823      Ashe
3  3 -80.97732  35.97633     Surry
4  4 -80.80150  36.20300 Currituck
5  5 -81.26429  36.23899      Ashe
6  6 -81.13721  36.31100     Surry

I think this might be possible by creating a "lookup/merge" script which takes each randomly simulated point and sees which location this point falls within?

Option 2: Or, perhaps there might be a more direct way to do this. For example, suppose from the shapefile, I could find out every pair of points that was located in "Ashe" - I could then just randomly sample these points and directly make a list of points within Ashe.

Can someone please help me in doing this?

Thank you!

Upvotes: 1

Views: 175

Answers (1)

nniloc
nniloc

Reputation: 4243

One option is to use the function sf::st_sample.

library(dplyr)
library(sf)

nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>% 
  st_transform(st_crs(4326)) %>% 
  st_cast('POLYGON')
#> Warning in st_cast.sf(., "POLYGON"): repeating attributes for all sub-geometries
#> for which they may not be constant

# use st_sample on Ashe to generate 10 random points
pts <- sf::st_sample(nc[1, ], 10)

# plot it
plot(st_geometry(nc[1:3, ]))
plot(pts, add = T, col = 'red')

Created on 2022-11-16 by the reprex package (v2.0.0)

Upvotes: 1

Related Questions