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 ---#
(<- readRDS('Data/well_registration.rds')
wells )
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 ---#
<- st_as_sf(wells, coords = c("longdd","latdd"))
wells_sf
#--- 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
bbox: 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 know48 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 ---#
<- st_set_crs(wells_sf, 4269)
wells_sf
#--- 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
bbox: xmin: -102.6249 ymin: 40.69824 xmax: -99.58401 ymax: 41.11699
geographic 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)
Yes, YOU need to know the CRS of your data.↩︎