2.5 Units
An important feature of CRSs is that they contain information about spatial units.Clearly, it is vital to know whether a house’s measurements are in feet or meters, and the same applies to maps.It is good cartographic practice to add a scale bar onto maps to demonstrate the relationship between distances on the page or screen and distances on the ground.Likewise, it is important to formally specify the units in which the geometry data or pixels are measured to provide context, and ensure that subsequent calculations are done in context.
A novel feature of geometry data in sf
objects is that they have native support for units.This means that distance, area and other geometric calculations in sf return values that come with a units
attribute, defined by the units package (Pebesma, Mailund, and Hiebert 2016).This is advantageous, preventing confusion caused by different units (most CRSs use meters, some use feet) and providing information on dimensionality.This is demonstrated in the code chunk below, which calculates the area of Luxembourg:
luxembourg = world[world$name_long == "Luxembourg", ]
st_area(luxembourg)
#> 2.41e+09 [m^2]
The output is in units of square meters (m2), showing that the result represents two-dimensional space.This information, stored as an attribute (which interested readers can discover with attributes(st_area(luxembourg))
), can feed into subsequent calculations that use units, such as population density (which is measured in people per unit area, typically per km2).Reporting units prevents confusion.To take the Luxembourg example, if the units remained unspecified, one could incorrectly assume that the units were in hectares.To translate the huge number into a more digestible size, it is tempting to divide the results by a million (the number of square meters in a square kilometer):
st_area(luxembourg) / 1000000
#> 2414 [m^2]
However, the result is incorrectly given again as square meters.The solution is to set the correct units with the units package:
units::set_units(st_area(luxembourg), km^2)
#> 2414 [km^2]
Units are of equal importance in the case of raster data.However, so far sf is the only spatial package that supports units, meaning that people working on raster data should approach changes in the units of analysis (for example, converting pixel widths from imperial to decimal units) with care.The new_raster
object (see above) uses a WGS84 projection with decimal degrees as units.Consequently, its resolution is also given in decimal degrees but you have to know it, since the res()
function simply returns a numeric vector.
res(new_raster)
#> [1] 0.000833 0.000833
If we used the UTM projection, the units would change.
repr = projectRaster(new_raster, crs = "+init=epsg:26912")
res(repr)
#> [1] 0.000833 0.000833
Again, the res()
command gives back a numeric vector without any unit, forcing us to know that the unit of the UTM projection is meters.