2.6 Turning a data.frame of points into an sf

Often times, you have a dataset with geographic coordinates as variables in a csv or other formats, which would not be recognized as a spatial dataset by R immediately when it is read into R. In this case, you need to identify which variables represent the geographic coordinates from the data set, and create an sf yourself. Fortunately, it is easy to do so using the st_as_sf() function. Let’s first read a dataset (irrigation wells in Nebraska):

#--- read irrigation well registration data ---#
(
  wells <- readRDS("Data/well_registration.rds")
)

        wellid ownerid        nrdname acres  regdate section     longdd
     1:      2  106106 Central Platte   160 12/30/55       3  -99.58401
     2:      3   14133   South Platte    46  4/29/31       8 -102.62495
     3:      4   14133   South Platte    46  4/29/31       8 -102.62495
     4:      5   14133   South Platte    46  4/29/31       8 -102.62495
     5:      6   15837 Central Platte   160  8/29/32      20  -99.62580
    ---                                                                
105818: 244568  135045 Upper Big Blue    NA  8/26/16      30  -97.58872
105819: 244569  105428    Little Blue    NA  8/26/16      24  -97.60752
105820: 244570  135045 Upper Big Blue    NA  8/26/16      30  -97.58294
105821: 244571  135045 Upper Big Blue    NA  8/26/16      25  -97.59775
105822: 244572  105428    Little Blue    NA  8/26/16      15  -97.64086
           latdd
     1: 40.69825
     2: 41.11699
     3: 41.11699
     4: 41.11699
     5: 40.73268
    ---         
105818: 40.89017
105819: 40.13257
105820: 40.88722
105821: 40.89639
105822: 40.13380

#--- check the class ---#
class(wells)

[1] "data.table" "data.frame"

As you can see the data is not an sf object. In this dataset, longdd and latdd represent longitude and latitude, respectively. We now turn the dataset into an sf object:

#--- recognize it as an sf ---#
wells_sf <- st_as_sf(wells, coords = c("longdd", "latdd"))

#--- take a look at the data ---#
head(wells_sf[, 1:5])

Simple feature collection with 6 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -102.6249 ymin: 40.69824 xmax: -99.58401 ymax: 41.11699
CRS:           NA
  wellid ownerid        nrdname acres  regdate                   geometry
1      2  106106 Central Platte   160 12/30/55 POINT (-99.58401 40.69825)
2      3   14133   South Platte    46  4/29/31 POINT (-102.6249 41.11699)
3      4   14133   South Platte    46  4/29/31 POINT (-102.6249 41.11699)
4      5   14133   South Platte    46  4/29/31 POINT (-102.6249 41.11699)
5      6   15837 Central Platte   160  8/29/32  POINT (-99.6258 40.73268)
6      7   90248 Central Platte   120  2/15/35 POINT (-99.64524 40.73164)

Note that the CRS of wells_sf is NA. Obviously, \(R\) does not know the reference system without you telling it. We know⁴⁸ that the geographic coordinates in the wells data is NAD 83 (\(epsg=4269\)) for this dataset. So, we can assign the right CRS using either st_set_crs() or st_crs().

#--- set CRS ---#
wells_sf <- st_set_crs(wells_sf, 4269)

#--- or this ---#
st_crs(wells_sf) <- 4269

#--- see the change ---#
head(wells_sf[, 1:5])

Simple feature collection with 6 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -102.6249 ymin: 40.69824 xmax: -99.58401 ymax: 41.11699
Geodetic CRS:  NAD83
  wellid ownerid        nrdname acres  regdate                   geometry
1      2  106106 Central Platte   160 12/30/55 POINT (-99.58401 40.69825)
2      3   14133   South Platte    46  4/29/31 POINT (-102.6249 41.11699)
3      4   14133   South Platte    46  4/29/31 POINT (-102.6249 41.11699)
4      5   14133   South Platte    46  4/29/31 POINT (-102.6249 41.11699)
5      6   15837 Central Platte   160  8/29/32  POINT (-99.6258 40.73268)
6      7   90248 Central Platte   120  2/15/35 POINT (-99.64524 40.73164)

---

48. Yes, YOU need to know the CRS of your data.↩︎