Welcome to the spatial half of the course. The objective is narrow and practical: to use R as your GIS, rather than clicking around in ArcGIS. This first deck is about vector data, meaning points, lines and polygons, handled through the sf package. Have a look at the table of contents, because it is a long deck and it moves from the inside out. We start with what an sf object actually is, then how to build one, then how to read and write them, and only at the end what you can compute with them. The prerequisite links on the right are worth a look if projections are unfamiliar. The full route is on the left. After the object structure and construction sections, you will cover file input and output, projection, quick plots, converting coordinate columns into sf, conversion to and from the older sp classes, ordinary data transformation, sticky geometry, geometric operations, and the failure modes you need for debugging. On the right, the prerequisite links refresh ggplot two and dplyr, which matter because sf plotting and manipulation build directly on them. The related links point beyond today’s basics to publication maps, topological subsetting, spatial joins, and vector-raster tasks such as cropping, masking, and extracting values. Treat those as the next steps after you can confidently inspect and manipulate the vector objects in this deck.
The objectives of this chapter are to learn how to use R as GIS, specifically how to handle vector spatial data.
sfsfg)sfc) and simple feature (sf)sfsp objectssfggplot2Two practical notes before we start. The navigation aids are the usual ones: stacked lines bottom left for a table of contents, letter o for a panel view. But the second half of this slide matters more than usual here. Spatial packages are large, and the first time you load one of these decks your browser has to download sf and its dependencies before anything will run. That takes a while, and there is no progress bar worth watching. So open the deck, give it a minute, and judge readiness by whether a cell actually produces output rather than by any status message. The pale blue area is the editable code area. Run Code evaluates the whole cell, while selecting only the lines you want and pressing Command plus Enter on a Mac, or Control plus Enter on Windows, evaluates just that selection. The two-sheet icon copies the cell so you can paste it into R on your own computer. The reload icon immediately to its left restores the original example after you experiment. Those controls are there so you can change an argument, run the result, and still return to the lecture version without reloading the entire deck. One more, at the right end of the toolbar: the eye icon hides that cell’s output and a second click brings it back. Some of these results run long, and once you have read one it is just pushing the rest of the slide out of view.
Click on the three horizontally stacked lines at the bottom left corner of the slide, then you will see table of contents, and you can jump to the section you want
Hit letter “o” on your keyboard and you will have a panel view of all the slides
Two things I am assuming, and both links are worth following if you are unsure. First, that you know roughly what a coordinate reference system is, and what projection means: the business of representing a curved earth on a flat plane, and the fact that you cannot do it without distorting something. Second, that you know the difference between vector and raster. Vector is points, lines and polygons with coordinates. Raster is a grid of cells. This deck is entirely vector; the next one is raster. If the CRS material is new, read that first link before going further, because projections come up repeatedly.
What the sf package gives you, and the first bullet is the important idea. A simple feature holds the geographic information and the attributes together in a data-frame-compatible structure, which is why sf feels natural if you already know data frames. The older sp package also stored attributes and geometry together, but used separate S4 spatial data-frame classes. Then the list underneath is a preview of the whole spatial section of this course: projection, reading and writing, non-interactive operations like buffers and areas, and interactive ones like spatial joins. Almost everything you will need for applied work is in that list. The remaining examples on screen make those categories concrete. Non-interactive operations change or measure one object, including buffers, areas, and distances. Interactive operations use spatial relationships between objects, such as subsetting by location or extracting values from features that intersect. Keeping all of these behind one sf-compatible structure is why you can move from ordinary attribute work to geographic work without rebuilding the dataset for every task.
The sf package provides a simple way of storing geographic information and the attributes of the geographic units in a single dataset called simple feature (sf).
The sf package allows you to do almost all the spatial operations you would need for your research
The dataset for most of this deck. It is North Carolina county boundaries, and it ships inside the sf package itself, which is why the path goes through system dot file rather than pointing at your disk. That makes it a convenient teaching dataset: everybody who has sf has this. We keep three columns, area, name and FIPS code, so the printout stays readable. Then run the class check on the second cell and look carefully at what comes back. It is not one class, it is two, and the next slide is about why. System dot file receives shape slash nc dot s h p and package equals sf, so it returns the installed path to that bundled shapefile. St read imports the path, and percent-greater-than-percent passes the resulting sf to dplyr double-colon select. The namespace prefix identifies dplyr explicitly, while AREA, NAME, and FIPS are the three attributes retained alongside the sticky geometry. The left arrow stores the result as nc. Class of nc then reports sf and data frame, reflecting the spatial layer built on top of a familiar tabular object.
Read the North Carolina county boundary data:
Check the class:
sfHere is the object itself, and the annotations on the right are worth reading line by line against the printout. The header tells you there are a hundred features and three fields. Then the body looks exactly like a data frame, rows as counties, columns as attributes, with one exception: that last column called geometry. That is the whole design. An sf object is a data frame that happens to carry a geometry column, which is why every dplyr verb you already know still works on it. Look at the Ashe County example in the callout and match it against the first row.
The first line tells you this is a simple feature (sf) object with 100 features and 3 attributes (fields)
So, an sf object looks just like a data.frame where rows represent observation units and columns represent attributes, except for a special column named geometry
The geometry column stores the geographic information of the observation units (here, county)
Example
Ashe County (1st row) has area of 0.114, FIPS code of 37009, and so on. And the entry in geometry column at the first row represents the geographic information of Ashe County.
Now zoom in on the geometry column on its own. Each entry is a simple feature geometry, an sfg, holding the shape of one county. Read the bullets on the right for the vocabulary, because these three abbreviations recur throughout the deck and they nest inside each other. An sfg is one shape. A collection of sfgs forms a simple feature geometry column, an sfc. And an sfc attached to a data frame is what makes an sf. Note also that these are all MULTIPOLYGON, because a county can consist of several disconnected pieces. The code uses dplyr double-colon select to ask nc for geometry alone. The namespace prefix says exactly which select function to use. Sticky geometry means this remains an sf view, but the one-column printout makes the hierarchy visible: each row contains one sfg, the full column is the sfc, and that sfc is the spatial component of nc.
An element of the geometry column is a simple feature geometry (sfg).
In general, sfg represents the geographic information of a single geometric feature (here, county).
There are different types of sfgs (POINT, LINESTRING, POLYGON, MULTIPOLYGON, etc)
In this example, all the sfgs are of type MULTIPOLYGON
A collection of multiple sfgs as a column is called simple feature geometry column (sfc), which can make a geometry column in an sf object
And one level deeper still, to see that there is nothing mysterious at the bottom. All those brackets are peeling the structure apart until we reach the raw numbers, and what we find is a matrix. Each row is a point, the first column is longitude, the second is latitude, and connecting the points in order traces the county boundary. That is all a polygon is. It is worth seeing this once, because it demystifies everything that follows: spatial operations are arithmetic on coordinates, not magic, and when something goes wrong it usually went wrong at this level. Read the extraction from left to right. Bracket one comma keeps the first county as an sf row. St geometry takes its geometry column, and the three successive double-bracket-one operations descend through the first feature, the first polygon, and the first ring of that multipolygon. Head with ten then limits the printed coordinate matrix to its first ten rows so you can inspect it without flooding the slide.
Let’s see what an sfg is made of.
Each row represents a point
Points are stored in a matrix format
Connecting all the points forms a polygon
sfg)The vocabulary of geometry types, and the examples in brackets are the useful part. A POINT has no area, so it suits a well, a city treated as a location, or a sampling site. A LINESTRING is a path, like a river reach. POLYGON is anything with area, like a county. Then the MULTI versions exist for a single thing made of several pieces: a river with tributaries, or a country with islands. That last one is why the North Carolina data came back as MULTIPOLYGON rather than POLYGON. Some of those counties include islands off the coast.
Some of the most common types of spatial objects represented by sfg are the following:
POINT: area-less feature that represents a point (e.g., well, city, farmland)
LINESTRING: (e.g., a tributary of a river)
MULTILINESTRING: (e.g., river with more than one tributary)
POLYGON: geometry with a positive area (e.g., county, state, country)
MULTIPOLYGON: collection of polygons to represent a single object (e.g., countries with islands: U.S., Japan, etc)
Starting with the simplest, and building one by hand is the point of the next few tabs. A POINT is just two numbers, so st point of a vector of two values makes one. Run it and look at how it prints: sf shows you the type in capitals followed by the coordinates. Then run the class check. The result has dimension class XY, geometry type POINT, and base geometry class sfg. Every geometry you make is an sfg of some particular type, and functions that work on sfg objects work on all of them regardless of type. The inner c combines two and one into the x-y coordinate vector, and the outer parentheses make R print the point at the same time that the left arrow assigns it to a point. The code-track value of zero point five gives the editor and its result equal shares of the live cell. That is only a display choice for these notes; it does not change the geometry R creates.
POINT is the simplest geometry type and is represented by a vector of two numeric values. An example below shows how a POINT feature can be made from scratch:
The st_point() function creates a POINT object when supplied with a vector of two numeric values. If you check the class of the newly created object,
you can see that it’s indeed a POINT object. But, it’s also an sfg object. So, a_point is an sfg object of type POINT.
A line is a sequence of points, so it needs a matrix rather than a vector. The first cell builds that matrix with rbind, four points, each a pair of coordinates. Then st linestring turns it into a geometry. Run the plot at the bottom and compare it against the numbers in the matrix. Each consecutive pair of points is joined by a straight line, in the order they appear. That ordering matters: shuffle the rows of the matrix and you get a different, probably self-crossing, shape from exactly the same set of points. The c calls make the four coordinate pairs, and rbind stacks those pairs as rows. The sf double-colon prefix names the package explicitly, so the function can be found even if sf has not been attached. Class confirms that the result is an X-Y LINESTRING and an sfg, while plot dispatches to the geometry’s own plotting method. Each code-track value of zero point five divides the live editor and output evenly on screen without affecting the R calculation.
A LINESTRING object is represented by a sequence of points in a matrix:
You can turn the matrix into a LINESTRING using sf::st_linestring():
Let’s plot it.
As you can see, each pair of consecutive points in the matrix are connected by a straight line to form a line.
A polygon is also a matrix of points, with one extra requirement stated in the sentence above the code: the first and last points must be identical, so that the boundary closes. Look at the matrix and check that, the first row and the last row are both zero zero. If you leave that out, you have a line rather than a region. Then note that st polygon takes the matrix wrapped in a list, which looks like unnecessary ceremony until the next tab, where the list is exactly what lets a polygon have more than one ring. The six coordinate pairs trace the exterior in row order, and class confirms that wrapping p one in the list produced a POLYGON sfg. In the last cell, plot draws that geometry and col equals red sets the polygon’s display colour. Code-track zero point four gives the editor forty percent of the side-by-side cell and leaves the larger share for its output. That layout option changes only the notes, not the object.
Just like the LINESTRING object we created earlier, a POLYGON is represented by a collection of points.
However, the first and last points in the matrix have to be the same to form a polygon
You can turn the matrix into a POLYGON using st_polygon(), which takes a matrix in a list() :
Let’s plot it.
And here is why the list mattered. A polygon can have holes, and the rule is in the first sentence: the first matrix in the list is the outer boundary, and every matrix after it is a hole cut out of it. So we reuse p one as the exterior and add p two, a small triangle inside it, as a hole. Run the plot and you can see the gap. This is not an exotic case. Any polygon with a lake in it, or a country surrounding an enclave, needs exactly this structure, and it is the reason st polygon has the interface it does. P two repeats one-one as its first and last coordinate so the triangular hole is itself a closed ring. List of p one comma p two preserves the order that gives the rings their roles. Plot then fills the valid polygon area red while leaving that interior ring empty. Code-track zero point four reserves forty percent of each live cell for the code and the rest for the printed or graphical result.
A POLYGON can have holes in it. The first matrix of a list becomes the exterior ring, and all the subsequent matrices will be holes within the exterior ring.
Let’s plot it.
One more level of nesting, for a single feature made of several separate polygons. The structure is a list of lists: each inner list is one polygon, in the same first-ring-then-holes format from the previous tab. Read the example carefully. The first inner list is the polygon with its hole, the second is a brand new polygon off to the side, and together they form one MULTIPOLYGON. Run the plot and you will see two disconnected shapes that R nonetheless treats as a single geometry. That is precisely how a county with an offshore island is stored. Rbind creates p three from five coordinate pairs, with four-zero repeated to close this second polygon. Autorun true makes that preparatory cell execute when the tab loads, so p three exists before the next cell uses it. The surrounding parentheses print p three while assigning it, and the next surrounding parentheses do the same for the completed multipolygon. Finally, plot draws both pieces in red. Code-track zero point five or zero point four controls how much horizontal room the editor receives in those live cells and has no effect on the list nesting.
To create a MULTIPOLYGON object you create a list of lists of matrices, with each inner list representing a polygon.
You supply a list of lists of matrices to the st_multipolygon() function to make a MULTIPOLYGON object.
Each of list(p1,p2) and list(p3) represents a polygon.
Your turn, and these three are deliberately unguided. Build a POINT, a LINESTRING and a POLYGON from scratch. You have seen all three, so this is about whether the structural rules have stuck rather than whether you can recall function names. Two things to watch. A line needs a matrix, not a vector, even if you only want two points. And a polygon needs its first and last rows to match, or sf will object. Plot each one as you make it, because seeing the shape is the quickest way to catch a coordinate you typed wrong.
POINTLINESTRINGPOLYGONsfc) and simple feature (sf)Now upward from single geometries to a column of them. An sfg is one shape; an sfc is a list-column of shapes. st sfc is what builds one, and you hand it a list of sfg objects. Look at what we are combining here, though, because it is unusual: a point, a line, a polygon and a multipolygon all in one sfc. Real data almost never mixes types like this, but it demonstrates that an sfc is genuinely just a list and does not require its contents to be uniform. Run the class check and note that a mixed geometry column has classes sfc GEOMETRY and sfc. Inspect individual elements to see their specific geometry types. The sf double-colon prefix makes the source package explicit. The left arrow stores the new column as sfc ex, while the surrounding parentheses print it immediately so you can see all four entries. Class then reports the common container classes rather than replacing the specific class carried by each element. Code-track zero point five gives the code and result equal width in both live cells.
sfg is an object class that represents a single spatial object.
We can combine multiple sfgs as a list to create a simple feature geometry list-column (sfc).
To make a simple feature geometry list-column (sfc), you can simply supply a list of sfg to the st_sfc() function as follows:
Check its class:
And the final step to a full sf object, in two nested tabs. First you make an ordinary data frame of attributes and assign the sfc to a column of it. Run the class check at that point and note what it says: still just a data frame. R does not infer spatialness from a column that happens to contain geometries. Then the second tab does the registration with st as sf, and only after that does the object report itself as an sf. That two-step is worth seeing, because it tells you what st as sf actually does. It sets metadata, not geometry. In the first tab, data frame creates four name values, A through D, so there is one attribute row for each of the four geometries. Assigning sfc ex to the dollar-geometry column keeps that row-for-shape alignment. In the second tab, st as sf recognises that existing sfc column as the active geometry column, and the final class call shows both sf and data-frame inheritance. Each code-track value of zero point four leaves forty percent of the live cell for code. The parentheses around the assignment print sf ex as it is created, which lets you inspect registration and values in one run.
To create an sf object, you first add an sfc as a column to a data.frame.
At this point, it is not yet recognized as an sf by R.
You can register it as an sf object using st_as_sf().
As you can see sf_ex is now recognized also as an sf object.
Two exercises building on each other. The first asks you to make an sfc from the point and polygon you created earlier, which is just st sfc over a list of the two. The second is the more interesting one: turn that sfc into a proper sf object with an id column, giving the point id one and the polygon id two. So you need a data frame with that id, the sfc attached as a column, and then the registration step. If you find yourself unsure whether you have an sf yet, check the class. That is what the previous tab was teaching you. The two bare object names in Exercise two point one print your existing point and polygon so you can confirm what you are combining before you call st sfc. In Exercise two point two, keep their order aligned with the id vector: the point is the first geometry and receives one, while the polygon is second and receives two. After attaching the geometry column, st as sf performs the same registration you just saw. Print the result or call class to verify both the row alignment and the sf class.
Create an sfc using the POINT and POLYGON you made earlier.
Create an sf object using the sfc object you created in the previous exercise, where the additional variable in the sf object is id with the POINT and POLYGON assigned id = 1 and id = 2, respectively.
A section about file formats, and the first bullet explains why we have to care. Most of the GIS world runs on ArcGIS, whose native format is the shapefile, so your collaborators will send you shapefiles and much public data is published as shapefiles. Then read the list of components, because it is the practical thing to remember. A shapefile requires the dot s h p, dot s h x, and dot d b f files. The dot p r j file is optional, but it should normally accompany them to preserve the coordinate reference system. Sending somebody the dot s h p on its own gives them something unreadable, which is exactly the annoyance the rest of this section goes on to complain about.
The vast majority of people still use ArcGIS software to handle spatial data, which has its own system of storing spatial data called shapefile system.
A shapefile requires three files, and a fourth should normally accompany them:
.shp: stores geometry information (like sfg).shx: a positional index into the .shp.dbf: the attribute table (everything that is not geometry).prj: optional, but preserves CRS information Chances are that your collaborators use shapefiles.
There are many GIS data online that are available only as shapefiles.
So, it is important to learn how to read and write shapefiles
Reading one in is a single function, st read, pointed at the dot s h p file. Note that you name only that one file even though the format requires three; sf finds the dot s h x and dot d b f siblings itself, provided they are sitting in the same folder with the same base name. The optional dot p r j file should travel with them so the coordinate reference system is preserved. That is also the most common way this goes wrong, so if a read fails, check that all the pieces travelled together. Then do the exercise in the callout with the download link. Reading somebody else’s shapefile off your own disk is the thing you will actually do first in a real project. In the syntax tab, file path is the argument that locates the dot s h p file, and st read returns the imported data as an sf object. The example assigns that result to nc imported and uses a path relative to the working directory, so Data is a folder inside the current project rather than an absolute location on this computer. Eval false keeps both displayed snippets from running while the deck renders because the example files live on your computer, not inside this browser session. Download every supporting file in the callout, copy the dot s h p path, and substitute that path in the same call.
You can use sf::st_read() to read a shapefile. It reads in a shapefile and turns the data into an sf object.
file_path: the path to the shapefile.Here, a file named nc_practice.shp (along with its supporting files) is read from the Data folder.
Try yourself
sf::st_read()Writing is the mirror image, st write, taking the object and a path. The one argument worth knowing is append equals FALSE, which is explained underneath. Without it, writing over an existing shapefile fails rather than replacing it, which is a sensible default that becomes maddening when you are iterating on a script. Setting it to FALSE says overwrite. Then the exercise asks you to write back out the file you just read in, which is worth actually doing, because it is how you will discover whether you had all four component files in the first place. The first positional argument is the sf object to export and the second is the destination path. In the example, nc imported is written under the new base name nc exported inside the project’s Data folder, and st write creates the required companion files alongside the dot s h p. Eval false keeps the illustrative syntax and example from touching a filesystem during rendering. When you try it yourself, include append equals false if that destination already exists and you intend to replace it.
You can use the sf::st_write() function to write an sf object to shape files.
append = FALSE forces writing the data when the shape files with the same name already existsThis code will export an sf object called nc_imported as nc_exported.shp (along with other supporting files) in the “Data” folder relative to the working directory.
Try yourself
sf::st_write() using whatever name you likeNow the argument that you should mostly not be using shapefiles at all, in four nested tabs. The motivation tab is blunt about it: the format is dominant for historical reasons rather than technical ones, and having four files per dataset is a nuisance. Then the alternatives. GeoJSON is one file and readable as text. GeoPackage is one file and technically the strongest of the three. And rds, if the data is only ever going to be read by R, which is both the simplest and usually the most compact. Read all four, then default to whichever your collaborators can open. In the GeoJSON tab, st write takes nc as the object and d-s-n as the destination name, here Data slash nc exported dot geojson. The filename extension tells the spatial driver which format to create, so there is no separate layer argument in this single-file example. St read reverses that operation and assigns the GeoJSON back to nc. Both chunks have eval false because they demonstrate local file operations rather than writing from the rendered lecture. The GeoPackage tab uses the same pattern with a dot g-p-k-g extension. It writes nc imported to one GeoPackage file and st read can restore it as nc. GeoPackage is especially useful when you want one broadly supported spatial file without shapefile’s collection of sidecars. Again, d-s-n names the destination and eval false prevents the examples from running during rendering. The final tab is different because rds is R’s own object format rather than a GIS interchange format. Save R-D-S takes the object first and its dot r-d-s path second, preserving the complete sf object in one compressed file. Read R-D-S takes that path and reconstructs the object, which the assignment names nc. Use this when the next reader also uses R; use GeoJSON or GeoPackage when another GIS must open the result. The callout’s storage point is the practical reason rds is attractive for large R-only workflows. These rds chunks are also marked not to evaluate, so reviewing or rendering the slide does not create or require any local files.
If your collaborators are using ArcGIS and demanding that they need a shapefile for their work, sure you can write to a shapefile.
But, there is really no need to work with the shapefile system if you are not using ArcGIS.
Basically, we are using the file system just because ArcGIS is the pioneer of GIS software and many people are still using it, but not because it is the best format available to store spatial objects.
Indeed, there are some limitations to shape files (see here).
But, first and foremost, it is annoying to have many files for a single spatial object.
A format that is increasingly popular is GeoJSON.
Unlike the shapefile system, it produces only a single file with .geojson extension.
GeoJSON files can also be read into ArcGIS.
Write
To write an sf object to a GeoJSON file, you simply give the file path to the dsn option (note that you do not use the layer option unlike the shape files case).
Read
You can use the sf::st_read() function to read a GeoJSON file like below:
One of the alternative data formats that is considered superior to the shapefile system is GeoPackage, which overcomes various limitations associated with shapefile.
Unlike the shapefile system, it produces only a single file with .gpkg extension.
GeoPackage files can also be read into ArcGIS.
Write
To write an sf object to a GeoPackage file, you simply give the file path to the dsn option (note that you do not use the layer option unlike the shape files case).
Read
You can use the sf::st_read() function to read a GeoPackage file like below:
Or better yet, if your collaborator uses R (or if it is only you who is going to use the data), then just save the sf object as an .rds file using saveRDS(), which can be of course read using readRDS().
Save
Read
Note
The use of rds files can be particularly attractive when the dataset is large because rds files are often more storage-efficient on disk than shapefiles.
A short motivation tab, and the reason it is short is that the need is entirely practical. Spatial operations such as joins require matching coordinate reference systems. Plotting is different: geom sf can plot layers with different known coordinate reference systems by transforming them automatically. So most of the time you reproject not because you have chosen a projection on its merits, but because something else you are working with is already in one. That is worth saying plainly, because it means the skill you need is not choosing projections well, it is checking and matching them reliably, which is what the next few tabs cover.
Spatial operations such as joins require matching coordinate reference systems (CRSs). By contrast, geom_sf() can plot layers with different known CRSs by transforming them automatically.
Checking what you have, with st crs. Run it and read the output, because it is verbose and mostly ignorable. What you are looking for is the last line, the EPSG code, which here is four two six seven. The bullets explain the surrounding format, Well Known Text, and the EPSG system: a catalogue of numbers standing in for full CRS definitions, so you can write one number instead of a paragraph. Note the caveat at the bottom, though. Only common systems have EPSG codes, so occasionally you will meet a CRS with no number and have to pass the full definition.
In order to check the current CRS for an sf object, you can use the sf::st_crs() function.
sf uses the Well Known Text format to store the coordinate reference system (CRS), which is one of many many formats to store CRS information (See here)
ID["EPSG", 4267] means that the EPSG code for this CRS is 4267
When you transform an sf using a different CRS, you can use its EPSG number if the CRS has an EPSG number
Reprojecting, in three nested tabs, and the third is the one to read twice. The mechanics are simple: st transform, the object, and an EPSG code. The example moves North Carolina into UTM zone seventeen north, which is the zone the state actually sits in, and then shows the geometry column before and after so you can see the coordinates change from degrees to metres. Then the Caveat tab shows st set crs, which sounds similar and is dangerous. It relabels the data without moving it, so your features stay where they were while claiming to be somewhere else. The How tab also allows a full Well Known Text definition in place of the EPSG number, and its snippet is marked not to evaluate because the names are spoken placeholders rather than runnable objects. In the example, st transform takes nc and target code two-six-nine-one-seven, then the left arrow stores the transformed copy as nc UTM while leaving nc unchanged. St crs checks the new metadata. The two select calls print only geometry from the original and transformed objects, which makes the coordinate change easy to compare even though sf’s sticky geometry rule would keep that column anyway. In the caveat, code two-six-nine-one-four is intentionally passed to st set crs to demonstrate relabelling. Do not use that result as reprojected data.
You can use sf::st_transform() to apply a different projection method to an sf object.
Syntax
Let’s transform (reproject) the data using NAD83 / UTM zone 17N CRS. Its EPSG number is 26917.
Let’s confirm the change in CRS:
Let’s compare the geometry column before and after the transformation (projection):
There is a function that sets CRS, namely sf::st_set_crs().
This function literally sets the CRS, but does not transform geometry accordingly unlike sf::st_transform().
So, doing this is a terrible mistake and the resulting sf object is no longer where it should be.
A small but genuinely useful trick. When you are reprojecting only in order to match another dataset, you do not need to look up its EPSG code at all. Pass st crs of that other object directly to st transform, and the two are guaranteed to agree. The example transforms back to where we started by asking for the CRS of the original. Prefer this over hard-coding numbers whenever the target is another dataset rather than a deliberate choice, because it cannot go subtly wrong the way a mistyped four-digit code can. Read the nested call from the inside out. St crs of nc extracts the complete target definition, and st transform applies that definition to nc UTM. The left arrow names the result nc UTM back to original, while the surrounding parentheses print it immediately. The next cell calls st crs on that result so you can compare it with the original definition. Code-track zero point three gives the code thirty percent of each side-by-side live cell and leaves more room for the long printed output.
You often need to change the CRS of an sf object when you interact (e.g., spatial subsetting, joining, etc) it with another sf object.
In such a case, you can extract the CRS of the other sf object using st_crs() and use it for transformation.
So, you do not need to find the EPSG of the CRS of the sf object you are interacting it with.
Let’s confirm the transformation:
Two short exercises on projection, using the fairway grid dataset, which is now loaded when the deck opens. The first is just st crs, to see what you have. The second asks you to find the EPSG code for NAD eighty three and transform to it, so you have to look the number up rather than being given it. That is deliberate: looking up EPSG codes is a routine part of spatial work and worth doing once here. And note this dataset is in Nebraska, which is why the later exercises use a different UTM zone from the North Carolina examples. For review, the geographic NAD eighty-three code you should find is four-two-six-nine. Use st transform with fairway grid as the object and that code as the target, assigning the returned sf if you want to keep it. The first exercise’s st crs result tells you what you are transforming from; the second changes the coordinates as well as their CRS metadata.
Check the CRS of fairway_grid.
Find the EPSG code for NAD 83, and change the CRS of fairway_grid to NAD 83 using the EPSG code.
The fastest possible look at spatial data: plot, with nothing else. Run it and note what you get, which is not one map but several, one per attribute column, each shaded by that column’s values. That is occasionally what you want and often surprising. The value of this is purely diagnostic. When you have just read a file or built a geometry, plot tells you in one keystroke whether the shapes are where you expect. It is not for making anything anybody else will see, and the note points forward to the ggplot deck for that.
The easiest way to visualize an sf object is to use plot():
plot() creates a map for each variable where the spatial units are color-differentiated based on the values of the variable
We will learn how to create more elaborate maps that are of publication-quality using the ggplot2 package later
One more viewing option, and note the instruction above the code: run this one on your own computer rather than here. mapView gives you an interactive map, a real slippy map you can pan and zoom, with a basemap underneath and your features on top, and clicking a feature shows its attributes. That is enormously useful for checking whether your data is where you think it is, which plot cannot really tell you. It does not run in these slides because the package is not loaded here, so copy it across and try it in RStudio. The first line uses sf double-colon st read, so the package source is explicit, and system dot file asks the installed sf package for its bundled shape slash nc dot s h p path. Package equals sf tells system dot file where to search, while quiet equals true suppresses the routine driver report. The next line calls mapview double-colon map view so you can use that package function without attaching the package. Eval false is what keeps the cell from executing in the deck. Message false and warning false suppress package chatter if you run it in a compatible render, and out-width eighty percent controls the displayed map width. Those chunk options affect execution or layout, not the two R objects.
Sometimes it is useful to be able to tell where certain spatial objects are and what values are associated with them on a map.
The mapView() function from the mapview package can create an interactive map where you can point to a spatial object and the associated information is revealed on the map.
Run the following codes on your computer.
data.frame of points into an sfA very common situation, and worth recognising when you are in it. Somebody sends you a csv with latitude and longitude columns. R reads it as an ordinary data frame, because nothing about a numeric column says it is a coordinate. So the data is spatial in principle and not in practice, and no spatial function will work on it until you say which columns are the coordinates. That conversion is one function call, which the next few tabs cover, but you have to know to make it, and you have to know the CRS to supply.
Often times, you have a dataset with geographic coordinates as variables in a csv or other formats
It would not be recognized immediately as a spatial dataset by R 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 sf::st_as_sf() function.
The dataset for this section: registered irrigation wells in Nebraska. Run the class check and note that it is a plain data frame, not an sf. The coordinates are sitting in two ordinary columns called longdd and latdd. Then read the callout, because it states the thing nobody can work out for you. You have to know your data’s coordinate reference system. It is not recoverable from the numbers, it comes from whoever produced the file, and here it is NAD eighty three, EPSG four two six nine. Guessing this wrong puts your data in the wrong place entirely. Data of wells N-E loads the packaged object, class verifies its current data-frame class, and head prints the first six records so you can identify long-d-d and lat-d-d among the attributes. Those names stand for decimal-degree longitude and latitude. Each code-track value of zero point five gives the editable code and its output equal horizontal space. The callout says G-R-S slash C-R-S on screen; the operational information you need here is the known coordinate reference system and its EPSG code.
Let’s get a dataset (irrigation wells in Nebraska) to work with:
wells_ne is a data.frame and has longdd and latdd representing longitude and latitude, respectively. Note that it is NOT an sf object.
Important
And the conversion, with st as sf. Three things to supply: the data, the coordinate columns, and the CRS. The one detail worth burning in is the order inside coords: longitude first, then latitude. That is x then y, which is the opposite of how people say it out loud, and getting it backwards is the classic beginner error. It does not throw an error either. Your points simply end up somewhere absurd, often in the ocean, which is at least a recognisable symptom once you have seen it happen. In the syntax, the first argument can be a data frame, tibble, or data table. Coords receives a character vector naming the x and y columns, and crs receives the definition those numbers already use. In the example, wells N-E is the input, long-d-d is x, lat-d-d is y, and four-two-six-nine identifies NAD eighty-three. By default, st as sf consumes those two coordinate columns and replaces them with a point geometry column while retaining the other attributes. The result is assigned to wells N-E sf. The setup-context cell prepares that object for later slides, and the visible cell repeats the conversion inside parentheses so it both assigns and prints the new sf. Code-track zero point four gives the code forty percent of the visible live cell.
We can turn a dataset (e.g., data.frame, tibble, data.table) into an sf object using sf::st_as_sf().
Your turn with a different dataset. The mower sensor data has columns called LAT and LNG, and you need to turn it into an sf and assign WGS eighty four using its EPSG code. Two things to work out for yourself. Which order those two columns go in, remembering the longitude-first rule from the previous tab, and what the EPSG code for WGS eighty four actually is. That last one is worth memorising, because it is the most common CRS you will meet: it is what GPS devices and web maps use. Data of mower sensor loads the exercise tibble. For review, the coordinate order is L-N-G then L-A-T, and the WGS eighty-four EPSG code is four-three-two-six. Pass mower sensor to st as sf, set coords to that longitude-then-latitude name vector, and set crs to four-three-two-six. The returned object is the spatial version, so assign it a new name if you want to preserve the original tibble for comparison.
Using the LAT (latitude) and LNG (longitude) columns, turn the tibble into an sf, and then assign the CRS of WGS 84 using its EPSG code.
sp objectsA short historical detour that occasionally matters. The sp package came before sf, by the same author, and while sf has replaced it for new work, a number of established packages still only accept sp objects. Two are named here, one for spatial econometrics and one for geographically weighted regression, and if your research needs either you will need this conversion. So this is not a section to learn deeply, it is one to remember exists, so that when a package rejects your sf object you know what to do about it.
The sp package is the predecessor of the sf package (developed by the same person)
There are many (older) packages that only accept spatial objects defined by the sp package
spdep: spatial econometricsGWmodel: runs geographically-weighted regressionIn that case, it is good to know how to convert an sf object to an sp object, vice versa.
Going from sf to sp, and the syntax is slightly odd because it uses base R’s as function with the string Spatial rather than a purpose-named converter. Run it and look at the class you get back. It will be something like SpatialPointsDataFrame, and note that the class name encodes the geometry type, which is a real difference from sf. In sf, one class holds any geometry. In sp, points, lines and polygons are different classes, which is part of why sf is more pleasant to work with once your data has mixed geometry types.
You can convert an sf object to its sp counterpart by as(sf_object, "Spatial")
And back the other way, which is simpler: st as sf accepts an sp object directly and works out the rest. Run it and check the class returns to sf and data frame. Worth noting that this is the same function you used earlier to turn a data frame of coordinates into an sf. It is doing the same conceptual job in both cases, taking something that holds spatial information in some other form and registering it as a simple feature, which is why it has the name it does.
You can convert an sp object to its sf counterpart by sf::st_as_sf(sp_object).
sfA short but important framing tab. Because an sf object is a data frame with an extra column, most of what you already know about manipulating data frames applies. That is not a small convenience, it is the main reason sf displaced its predecessor. You do not need to learn a separate spatial dialect for filtering rows or computing a new variable. The dplyr verbs you spent Chapter three on work here, and the next tab demonstrates them. Some verbs have geometry-specific behavior, such as sticky select and geometry-unioning summarize, which you need to understand.
An important feature of an sf object is that it is basically a data.frame with geometric information stored as a variable (column).
Most data-frame verbs work on an sf object.
Some verbs have geometry-specific behavior, such as sticky select() and geometry-unioning summarize().
dplyr verbs work well with sf
And the demonstration. Select, filter and mutate use the familiar data-frame syntax on an sf object. The one behaviour worth noticing is called out twice: the geometry column survives a select even when you did not ask for it. That is deliberate and it is what you want almost always, since a spatial object without its geometry is not spatial. It does mean select cannot be used to drop geometry, though. That takes a specific function, and the note at the bottom points out that none of this worked with sp objects. The first call selects well-id from wells N-E sf and demonstrates that geometry remains beside it. In the longer pipeline, each percent-greater-than-percent passes the current result to the next verb. Select keeps well-id plus the sticky geometry, filter keeps only rows whose well-id is greater than twenty thousand, and mutate replaces well-id with its old value plus twenty. The feature rows and their shapes remain aligned through every step. Writing dplyr double-colon before each verb makes its package source explicit, and the callout’s contrast with sp is why the data-frame compatibility of sf matters in practice.
The following code selects wellid variable using dplyr::select():
Notice that geometry column will be retained after dplyr::select() even if you did not tell R to keep it above.
Of course, you can apply other dplyr verbs just like you do with a data.frame. Here, let’s apply dplyr::select(), dplyr::filter(), and dplyr::mutate() in sequence using a piping operator.
Note
You cannot do this with the spatial objects defined by the sp package
Here is a property of sf objects that is mostly a blessing and occasionally a nuisance. The geometry column will not let go. Select two other columns and geometry comes along uninvited. Filter to ten rows and the right ten shapes follow. That is deliberate and it is the reason sf is safe to use: your shapes cannot drift out of alignment with your attributes, which is exactly the bug that plagued the older way of doing this. But it has consequences worth knowing, and the next few tabs are the ones that catch people. The escape hatch, when you want it, is st_drop_geometry.
An sf object is a data.frame with a geometry column that refuses to be left behind. Select other columns and it comes along. Filter rows and it follows. Group and summarise and it does something you may not have asked for.
That stickiness is the point of sf — your geometry can never drift out of sync with your attributes. But it has consequences, and the next three tabs cover the ones that actually catch people.
sf::st_drop_geometry() is the escape hatch. It returns an ordinary data.frame.
Group an sf and summarise it, and something happens that no other package would do: the geometries of each group are merged into a single shape. A hundred counties become two rows, and each row’s geometry is the union of the counties in its group. Look at the map, where the county lines inside each group have vanished. Whether that is wonderful or annoying depends entirely on what you wanted. If you were dissolving counties into regions, this is the tool and you did not even have to ask for it. If you only wanted a group mean, it is a lot of work you did not want. The first pipeline starts with nc. Mutate creates grp by repeating a and b to the number of rows in nc; length dot out makes the repetition stop at exactly that length, and the dot inside nrow refers to the piped nc object. Group by then forms the two groups. Summarize calculates mean area for each and, because the input is sf, unions every county geometry within the same group. The assignment stores those two dissolved features as nc grouped, and the final line prints them. Autorun true prepares that object as soon as the cell loads. The plot starts with an empty ggplot and adds one geom sf layer whose local data is nc grouped. Mapping grp to fill inside aes gives the two regions different fills and creates the legend. Alpha zero point five makes both fills half transparent, and theme void removes axes and background details so the dissolved boundaries are the focus. Code-track zero point four leaves more width for that map. The callout gives you the decision rule: keep geometry when you want the dissolve, and drop it before summarizing when you only want the grouped attribute result.
Group an sf and summarise it, and the geometries of each group are unioned into one shape:
Sometimes exactly what you want
One hundred counties became two rows, and each row’s geometry is the union of the counties in that group. If you wanted to dissolve counties into regions, this is a gift — it is how you do it, and you did not have to ask.
If you only wanted a mean, it is wasted work, and on a large object it is a lot of wasted work. Drop the geometry first and it becomes an ordinary dplyr operation.
And here is the price of that convenience. Run the block and compare the two timings — and be patient with the first one, because it is genuinely that slow. Collapsing rows means merging their shapes, so the version that keeps the geometry takes hundreds of times as long as the one that drops it first. Unioning polygons is real geometric computation; averaging a column is almost free. But be careful how you generalise from that, because it is easy to over-learn. Filter, mutate, select and arrange are not slow on an sf. They carry the geometry column along but they never recompute it, so they cost about what they would on a plain data frame. It is summarise specifically that hurts. So do not go dropping geometry everywhere out of superstition. Drop it before a summarise that does not need a shape, and join the result back afterwards on an ID column. System dot time measures the elapsed and processor time for each braced loop. In both loops, i runs from one through five so the same summary is repeated enough times to compare. The first pipeline groups nc by county name and computes m as the mean of area. Even though each name identifies one county here, sf’s summarize method still performs its geometry handling. The second pipeline inserts st drop geometry before the same group by and summarize, so the remainder operates on an ordinary data frame and never asks the geometry engine to union anything. The repeated result is not saved because only the timing comparison matters. Read the elapsed entry in the output for the wall-clock contrast, then apply the callout’s narrow lesson to summarize rather than blaming all dplyr verbs.
summarize() is where the stickiness costs you: collapsing rows means unioning their polygons, wanted or not.
It is summarize(), not dplyr
Keeping the geometry costs hundreds of times as long — about 400-500x in this browser, roughly 200x in RStudio — and the gap widens as the data grows.
filter(), mutate(), select() and arrange() never rebuild the geometry, they only carry it along, so they cost milliseconds either way. Drop geometry before a summarize() that does not need a shape, not defensively everywhere, and join the result back on an ID column.
This one is genuinely surprising the first time. When you join an sf to a plain data frame, the order of the arguments decides what class comes back. Put the sf first and you get an sf. Put the data frame first and you get a plain data frame, which still contains a column called geometry holding all the right shapes, so it prints perfectly happily and looks fine. But it is not an sf, so geom sf and many spatial operations require you to convert it first. The geometry column remains an sfc, which is why st as sf can convert it back. Put the sf first, or repair it afterwards with st as sf. The preparation cell builds county pop as a plain data frame. Its NAME column copies nc dollar NAME, and seq len of nrow nc creates the integers one through one hundred as a stand-in population value, one per county. Head prints the first six rows so you can inspect the join key and values. Code-track zero point five splits that cell evenly between code and output. Both left joins use by equals NAME, so NAME is the matching key and every row from the first argument is retained. In joined a, nc is first and county pop supplies pop, so sf’s method preserves the spatial class. In joined b, county pop is first and nc supplies the geometry-bearing rows, so ordinary data-frame dispatch determines the class. The outer parentheses print each assignment immediately, and code-track zero point five again gives each editor half the cell. The callout’s repair, st as sf of joined b, works because that second join changed the container class but did not destroy the sfc column.
When you join an sf to a plain data.frame, the order of the arguments decides what you get back.
sf first:
data.frame first:
The second one is not an sf
Put the data.frame first and you get a plain data.frame back. It still has a column called geometry holding the shapes, so it looks fine when you print it; but it is not an sf, and geom_sf() and many spatial operations require conversion.
Two ways out: put the sf first, which is the habit to build, or repair it afterwards with st_as_sf(joined_b), which works because the geometry column remains an sfc.
The last one, and it is the most disappointing to discover late. CSV has no way to represent a polygon, so if you write an sf to CSV the geometry is flattened into a text dump of the raw coordinate list. One county turns into a five-hundred-character line, and nothing on earth will read the shape back out of it. The rule is straightforward. If you want the attributes, drop the geometry and write that. If you want the spatial object, use st_write with a geopackage or geojson extension, which we covered in the reading and writing section earlier in this deck. The diagnostic pipeline takes row one of nc, converts that one-row sf to a plain data frame without dropping its geometry column, formats every value as text, and then uses substr from character one through ninety so the slide shows only the beginning of the enormous representation. That preview explains the literal list-of-coordinate text in the callout. St drop geometry is different from as data frame here: it deliberately removes the geometry before CSV export rather than merely removing the sf class.
Do not write.csv() an sf
CSV has no way to hold a polygon, so the geometry is flattened into text like
0.114,"Ashe","37009",list(list(c(-81.4727554321289, -81.5408401489258, ...One county becomes a 524-character line, and nothing can read the shape back. If you want the attributes, st_drop_geometry() first and write that. If you want the spatial object, use st_write() with a .gpkg or .geojson extension, as in the reading and writing section.
Buffers, in three nested tabs, and the pair of examples together make a point worth more than either alone. The How tab gives the syntax and the crucial hint: distance is in the units of your CRS, and you can check those with st crs of your object dollar units. Then the non-projected example buffers unprojected data. Look at the result and you will see a zig-zag edge, which the note explains and links out about. Then the projected example does it properly, transforming into the local UTM zone first. That is the recommendation to take away: project, then measure. The syntax says st buffer can expand points, lines, or polygons. Its first argument is the spatial object, and dist supplies the radius in the object’s working units. The syntax chunk is not evaluated because sf and distance are placeholders. In the first example, bracket one comma selects the first North Carolina county, dist equals two thousand requests the buffer, the left arrow names it nc buffer, and the surrounding parentheses print it. Warning true leaves any warning visible for this example, while code-track zero point five divides code and output evenly. The first ggplot draws the original county outline red, then draws nc buffer with a blue fill at alpha zero point three so you can still see the original through it. Theme void removes non-map decoration. In the projected version, st transform first converts that county to WGS eighty-four UTM zone seventeen north using code three-two-six-one-seven. The percent-greater-than-percent then passes the projected feature into st buffer, where two thousand now plainly means two thousand metres. Its plot again shows the original in red and the buffer in translucent blue. Geom sf can transform the two known CRSs for display, even though spatial analysis between mismatched CRSs would fail. The code-track options allocate editor width only, and the handoff from the zig-zag example to the projected tab is the workflow to keep: transform first, then buffer.
You can use sf::st_buffer() to create buffers of the specified length around points, lines, and polygons.
Syntax
dist: provide the distance in the unit of the CRS (run st_crs(sf)$units to get the unit)Let’s create a buffer of 2000 meter.
Here is what it looks like:
Yes, you see zig-zag. You can read up on this here. For now, I would recommend that you project first and then create a buffer.
Let’s first project and then create a buffer:
Here is what it looks like:
Computing areas, and the second nested tab is the one that catches people. st area gives you the area of every polygon, which is straightforward. But run the class check and note what comes back is not a plain number, it is a units object, carrying square metres around with it. Units objects support arithmetic and automatically convert compatible units, while rejecting dimensionally incompatible operations. The Caveat tab shows how to wrap the result in as numeric when a downstream operation specifically requires plain numeric values. Do that deliberately rather than by accident, because you are discarding the unit information. In the Example tab, mutate adds area by calling st area on all geometries in nc. The pipe then sends the result to select area. Because geometry is sticky, the printed nc still contains the geometry beside that one requested attribute. The assignment deliberately replaces nc with this narrower sf, and the surrounding parentheses print it. Code-track zero point five gives code and output equal space. In Caveat, class of nc dollar area exposes the units class. The next mutate recomputes the areas and pipes them through as numeric, which strips the attached square-metre metadata while keeping the magnitudes. Select again retains area and sticky geometry, and class now shows an ordinary numeric vector. Code-track zero point four leaves more room for those outputs. Keep the units version unless a later function actually demands numeric, because once you strip the class R can no longer protect you from mixing square metres with some other area unit.
You can use st_area() to calculate the area of all the polygons in an sf object.
By default, area calculated by st_area() is units.
Units objects support arithmetic with automatic conversion of compatible units. Drop the units like this only when a downstream operation specifically requires plain numeric values:
Centroids, in two nested tabs. st centroid replaces each polygon with a point at its center of mass, which can be useful for computing rough distances between regions. A centroid can fall outside a concave polygon, so use st point on surface when a label point must lie inside. Then the second tab pulls the raw coordinates out with st coordinates, giving you a plain matrix of x and y, and shows how to bind that onto the original object. That is the standard route when you need coordinates as ordinary columns, for example to pass them to something that is not spatially aware. The first call applies st centroid to every feature in nc, assigns the returned point sf to nc centroids, and prints it because the assignment is wrapped in parentheses. The plot then starts an empty ggplot, draws the original county geometries, and adds the centroid points as a second geom sf layer. That overlay lets you judge each point against its source polygon. Each code-track value of zero point four gives the editor forty percent of the cell. In Matrix of coordinates, the pipe applies st centroid once more to nc centroids, which are already points, then st coordinates converts those point geometries to an X-Y matrix, and head limits the display to six rows. The final line binds the coordinate matrix column-wise onto nc and pipes the combined object to head. Cbind is useful here because both objects have one row per feature in the same order. It does not perform a key-based match, so that row alignment is the condition that makes the result meaningful.
You can use st_centroid() to find the centroid of each of all the polygons in an sf object.
As you can see, st_centroid() returns an sf of centroids as points. The centroids look like this:
If you want longitude (X) and latitude (Y) of the centroids, you can further apply st_coordinates().
Of course, you can easily add the XY matrix to the original sf file using cbind():
Distances, and there are two things to notice. First the shape of the result: st distance between two sets of features gives you a full matrix, every feature in the first against every feature in the second, not a paired-up vector. Read the description above the code carefully so you know which index is which. Second, and this genuinely surprises people, read the callout. Even though this data is unprojected, in degrees, the distances come back in metres. sf computes them on the sphere rather than treating degrees as if they were a flat grid. The How tab’s first and second arguments are the two feature sets, and the displayed call is not evaluated because sf one and sf two are placeholders. In the example, rows one through five of nc are reduced to five centroids for the matrix rows, while rows six through fifteen become ten centroids for its columns. The result is therefore a five-by-ten units matrix, and entry i comma j measures the corresponding row centroid against the corresponding column centroid. Code-track zero point three leaves most of the live cell width for that matrix. The metre result is the useful longitude-latitude exception that the callout asks you to notice.
Syntax
This finds the distance between each of the points in sf_1 and each of the points in sf_2.
Get a matrix of distances whose [i,j] element is the distance between the ith sfg of st_centroid(nc[1:5, ]) and jth sfg of st_centroid(nc[6:15, ]):
Note
Notice that, even though nc is unprojected, distances returned are in meter.
The last operation: dissolving many geometries into one with st union. The example collapses a hundred counties into a single outline of the state, which is what you want when you need a boundary rather than its subdivisions. Note the warning in the text, though, because it is easy to trip over. What comes back is an sfc, not an sf. The attributes are gone, and reasonably so, since there is no single sensible answer for what the name or area column of a merged object should be. If you need attributes, you have to decide how to aggregate them yourself. The left arrow stores the unioned geometry as nc one, and class verifies the sfc result before plot draws its single dissolved outline. Because st union receives nc without a second object, it combines all geometries within that object rather than testing a relationship between two datasets. Each code-track value of zero point three gives the code a smaller share so the printed class and map have more room. A common next step, such as finding one centroid for the whole state, would operate on nc one rather than on the hundred county shapes.
Sometimes you want to combine all the geometries in a single sf. For example, you may want to get the centroid of North Carolina using nc (of course, you can alternatively get sf of NC state boundary, instead of counties in this case).
You can use sf::st_union() to achieve this. Note that the returned object is sfc, not sf.
Here is what it looks like.
Everything so far has shown you operations that work. This section is about the three ways they usually do not, because you will meet all three within about a week of using real spatial data, and none of them announces itself clearly. Different coordinate systems, broken geometries, and distances in units you did not expect. The reason to spend a slide on each is that the error messages are genuinely unhelpful, and recognising one on sight is the difference between a one-line fix and an afternoon of guessing where you went wrong.
Three things go wrong far more often than anything else, and all three produce error messages that do not obviously say what is wrong. Knowing them by sight turns an afternoon of confusion into a one-line fix.
This is far and away the most common failure, and it has a signature. Whenever two objects do not share a coordinate reference system, sf refuses to combine them and says st_crs of x equals st_crs of y is not TRUE. Learn that sentence by sight. It comes from intersection, from joins, from distance, from almost everything that takes two objects. The fix is always a transform, and the only judgement is which object to move. Move vector data rather than raster data, because re-projecting a raster resamples its values and you cannot undo it. The first live cell creates the mismatch deliberately. St transform makes nc UTM in WGS eighty-four UTM zone seventeen north with code three-two-six-one-seven, and comparing st crs of nc with st crs of nc UTM returns false. Autorun true prepares that object and result when the slide loads, while code-track zero point five splits editor and output evenly. The next cell passes the first five features from each object to st intersection. Those features describe the same place, but their coordinate numbers are expressed in different systems, so the operation stops before trying to intersect them. The callout’s fix transforms nc UTM back to the CRS extracted from nc inside the st intersection call. That snippet is marked not to evaluate because it is the repair pattern to read rather than another live failure. When both objects are vector data, choose the target CRS that makes sense for the analysis. The warning about rasters explains why the apparently symmetric choice is not symmetric when one input is a grid.
This is the most common failure by a wide margin, and it has a signature message:
Now try to combine them:
st_crs(x) == st_crs(y) is not TRUE
That exact wording is the one to memorise. sf refuses to combine two objects unless their CRS match, and you get the same message from st_intersection(), st_join(), st_distance() and most of their relatives.
The fix is always the same: transform one object to the other’s CRS.
Which one should move? Transform the vector data, not raster data, and when both are vector, transform to whichever CRS the analysis needs. Re-projecting a raster resamples its values and cannot be undone.
A geometry can be well-formed as data and still describe something impossible. The bowtie here crosses itself, and real downloaded shapefiles contain shapes like it routinely. Watch what happens when we ask for its area, because this is the important part: it returns zero, with no error and no warning, since the two triangles cancel. You would believe that number. Other operations do refuse, with messages about topology exceptions or self-intersections, and the wording varies by version. st_make_valid repairs it, but notice it returns two polygons where there was one. Repair changes your data. Start in A broken shape and read the construction from the inside out. Rbind stacks five coordinate pairs whose last row repeats zero-zero to close the ring. Their order traces two diagonals that cross at one-one. St polygon wraps that ring as a polygon, st sfc turns it into a geometry column and assigns projected CRS three-two-six-one-seven, and st sf attaches the name bowtie to make a one-row sf object. Autorun true builds it on load, and code-track zero point five gives code and output equal space. The ggplot layer uses bowtie as its data, fills it red, and sets alpha to zero point three so the crossing remains visible. Theme void removes axes and background. St is valid first returns a logical answer. With reason equals true, it returns the diagnostic text and the one-one crossing coordinate instead. That argument matters because false alone tells you only that something is broken, while the reason helps you locate what. Move to It answers anyway and st area returns zero units despite the visible shape. The callout explains why silence is dangerous here: validity must be checked rather than inferred from the absence of an error. In When it does complain, intersecting bowtie with itself forces a topology operation and may produce one of the listed messages. The exact words can vary, so recognise the family of topology and self-intersection complaints rather than memorising only one string. Finally, Repairing it assigns st make valid of bowtie to bowtie fixed. St is valid confirms the repair and st area gives the repaired area. The callout is part of the result, not a footnote: the repair interprets the crossing ring as two triangles, so the geometry becomes a multipolygon and its meaning changes. Inspect feature type, feature count, and area after any repair. For longitude-latitude data, follow the callout’s additional advice if the spherical repair remains invalid: project first, repair there, and inspect again.
A geometry can be syntactically fine and still describe something impossible — a boundary that crosses itself, for instance. Real downloaded shapefiles contain these routinely.
A “bowtie”: the outline crosses itself in the middle
Ask whether it is valid, and why not:
“Self-intersection[1 1]” — the coordinates where it crosses itself
Now watch carefully. Ask for its area:
It answered. The answer is wrong.
You get 0, with no error and no warning. The two triangles have equal area with opposite orientation, so they cancel exactly. An invalid geometry does not necessarily stop your code — it can quietly hand back a number you will believe.
Other operations do refuse:
What the complaint looks like
TopologyException, side location conflict, Self-intersection, Loop 0 is not valid, Edge N crosses edge M: the exact wording depends on which geometry engine handled the call, and differs between sf versions and between projected and unprojected data. They all mean that the geometry is invalid. Inspect the input geometry and preceding operations to determine where invalidity was introduced.
st_make_valid() repairs it:
Repairing changes the geometry
Look at what came back. The single self-crossing POLYGON is now a MULTIPOLYGON of two triangles, because that is the only sensible reading of the shape. So st_make_valid() is not a formality — it can change feature counts and areas, and you should look at what it did rather than piping straight through it.
One wrinkle worth knowing: on unprojected lon/lat data, st_make_valid() goes through the s2 spherical engine and does not always succeed. If it returns something still invalid, project the data first and repair it there.
The last one, and the sneakiest, because nothing errors at all. A number you hand to st_buffer is in the units of the coordinate system, and those are not always metres. Buffer this county by two thousand on a UTM projection and it grows to about fourteen hundred and thirty-one square kilometres. Buffer it by the same two thousand on the North Carolina state plane system, which is measured in survey feet, and the smaller six-hundred-and-ten-metre buffer grows it to about twelve hundred and twenty-five square kilometres. Convert the areas to common units before dropping their units, and check st crs units before passing any distance. Longitude and latitude, oddly, gives you metres. In Two CRS, two units, bracket one comma keeps the first county as nc one. St transform makes nc metres with UTM code three-two-six-one-seven and nc feet with state-plane code two-two-six-four. Dollar-units extracts the unit label from each st crs result, so you can verify metre versus U-S survey foot before measuring. Autorun true prepares both transformed objects, and code-track zero point five divides the live cell evenly. In Same number, different distance, both st buffer calls use dist equals two thousand, but each object interprets that number in its own CRS units. St area then returns a units-aware area for each buffer. Units double-colon set units converts each result to kilometres squared while it still knows the source unit, and only then does as numeric remove the units class. The outer c combines the two comparable values and names them metres and feet. Code-track zero point five again gives the editor and result equal widths. The important callout explains both the unequal areas and the safe order of operations. A raw square-foot value divided by one million is still not square kilometres. Convert first and strip second. Then read Lon slash lat is the exception as a deliberate contrast: for unprojected longitude-latitude objects, sf delegates to the s-two spherical engine and treats the buffer distance as metres. That convenience does not remove the habit of checking st crs of x dollar units for projected data.
A distance you pass to st_buffer() or st_is_within_distance() is expressed in the units of the CRS, and those are not always what you assume.
The same number now means two very different distances:
dist = 2000 is not a distance until you know the CRS
On the metre-based CRS the county grows to about 1,431 km². On the feet-based one it grows to about 1,225 km², because 2000 US survey feet is only about 610 metres. Convert units before dropping them; dividing a raw square-foot value by one million does not convert it to square kilometres.
Check st_crs(x)$units before passing any number as a distance. Plenty of US state plane systems are in feet, and nothing warns you.
Longitude/latitude is the helpful exception
You might expect a number to mean degrees on unprojected data. It does not. sf hands lon/lat work to the s2 spherical engine, which measures in metres, so st_buffer(nc, dist = 2000) on unprojected nc really does buffer by 2 km. That is convenient, and it is also why unprojected distances quietly work when you expected them to fail.
Now a sequence of four exercises that chain together, using the fairway grid dataset from a Nebraska golf course. Start here by loading it and looking at it, and do actually plot it before moving on, because each exercise builds on the one before and a mistake early will follow you through all four. The sequence mirrors the whole deck: project, then buffer, then find centroids, then measure distances between them. If you can do those four in order, you can do most of what applied vector work asks for.
Run the following code to get fairway_grid data. Then, inspect the data (e.g., by plotting) to get a good sense of what the data looks like.
The first exercise, in two parts. Plot the data, then transform it to NAD eighty three UTM zone fourteen north, whose code you are given. Note that fourteen north is correct here even though the North Carolina examples used seventeen north. This dataset is in Nebraska, and the whole point of UTM is that you use the zone your data is in. Getting that habit right is more valuable than remembering any particular number. Name the result exactly as the exercise asks, because the next three exercises all expect to find it under that name. The folded answer first calls plot on fairway grid for the diagnostic map. It then passes fairway grid and target code two-six-nine-one-four to st transform and assigns the returned sf to fairway grid UTM. Code-fold true keeps that solution collapsed until you choose to reveal it, and eval false prevents the answer cell from running for you. Those options preserve the exercise while still leaving a reviewable solution on the slide.
First, plot fairway_grid to get a sense of what the dataset looks like:
Transform fairway_grid so that its CRS is NAD 83/UTM zone 14N (its EPSG code is 26914) and name it fairway_grid_utm.
The second exercise: ten metre buffers around each grid polygon. This is where projecting in the previous exercise pays off, and it is worth pausing on why. Because the data is now in UTM, the units are metres, so distance ten means ten metres and nothing further is needed. Had you skipped the transform and buffered the unprojected data, you would have got the zig-zag edges from the buffers tab. So the order of these two exercises is itself the lesson, rather than an accident of how they happen to be written. The answer passes fairway grid UTM as the spatial input to st buffer, sets dist to ten, and assigns the returned buffer polygons to fairway grid buffers for the next exercise. Code-fold true hides the solution until requested, and eval false means revealing it does not execute it. Run your own live cell first, then compare the function arguments and object name.
Create buffers around the grid polygons in fairway_grid_utm where the radius of the buffer is 10 meter, and name it fairway_grid_buffers.
The third exercise: centroids of the buffers you just made. Mechanically this is one call to st centroid, so it should be quick. Worth thinking about what you expect before you run it, though. Buffering an irregular or asymmetric polygon can shift its centroid because the added area is not necessarily distributed symmetrically around the original centroid. Calculate and compare the centroids rather than assuming that buffering preserves their locations. If a shift looks unexpected, inspect the original shapes and the preceding operations. The answer calls st centroid on fairway grid buffers and assigns the point sf to buffers centroids, exactly the name the final exercise expects. Code-fold true makes that answer expandable, and eval false leaves execution to your live work area. Inspecting or plotting buffers centroids after the assignment is a useful check that there is one point for every buffered feature.
Find the centroid of each of the buffer polygons you created in Exercise 2, and then name it buffers_centroids.
And the last one: distances between those centroids. Note that this is st distance with a single argument rather than two, which computes every pairwise distance within one object. So what you get back is a square matrix with zeros down the diagonal, since every point is zero from itself. That diagonal is a useful sanity check. And because you projected back in exercise one, these distances are in metres from the CRS rather than computed on the sphere, which is exactly what you want for a dataset at the scale of a single field. The answer is simply st distance of buffers centroids. With n centroid features, its result has n rows and n columns, is symmetric across the diagonal, and carries distance units. Code-fold true keeps that one-line solution hidden until you open it, while eval false prevents the deck from calculating it automatically. Compare the matrix dimensions and units as well as the numeric entries when you check your work.
Calculate the distances between the centroids in buffers_centroids.