Geographic coordinate system vs projected coordinate system (Lat long vs UTM)
Both vector (point, line, polygon) and raster (pixels) data rely on x, y values to situate the data in 2D space to represent where it is on the earth’s surface. All coordinate systems establish a starting point on the earth (the origin) and measures from there.
There are different datums used to model the shape of the earth. It consists of a series of numbers that define the shape and size of the ellipsoid and it’s orientation in space. There are a large number of datums in use for various areas such as: The WGS84 datum, which is almost identical to the NAD83 datum used in North America and the ETRS89 datum used in Europe. You must select a datum that is best suited for your geographic location to give the best possible fit.
Map projections are mathematical equations used to flatten spherical objects. The projection formulas allow us to convert from lat/long coordinates (spherical system) to planar x, y coordinates (i.e. UTM). Lat/long coordinates are stored in degrees vs UTM the coordinates are stored in meters. Projected coordinate systems (PCS) are often better suited for large-scale applications because distances are often in meters. They are especially good for raster, because the pixel sizes are in meters rather than decimal degrees.
The Universal Transverse Mercator (UTM) projection uses a grid defined by Easting (vertical lines referenced by a zone’s central meridian) and Northings (horizontal lines measures from the equator) to locate coordinate pairs. The UTM system divides the world into 60 N-S zones, each representing 6° (e.g.,
360°/60 = zone width of 6° longitude). The zones are numbered from 1 to 60, starting at 180° W and ending at 180° E. In BC we are covered by UTM zone 8, 9, 10 and 11. The Okanagan occurs mostly in UTM zone 11, with a small portion in 10. Example of UTM coordinates: 11U 328058mE, 5534668mN.
When you are ‘projecting’ from a geographic to a projected coordinate systems (or vice versa) you are actually applying a mathematical transformation to the data to convert the coordinates from degrees to meters. This will ALWAYS distort the data slightly, which is why for a raster every time we project the data we also resample it slightly because the pixels are being shifted.
In R, things work a bit different. When you run: crs(raster), or any vector, it will return the proj4 string that mathematically defines how the coordinates are stored with the data. This is a proj4 string: “+proj=longlat +datum=WGS84 +no_defs“. It tells me my data are in a longlat (so GCS) projection that uses the World Geodetic Survey of 1984 survey (the standard for most GPS). For a raster in R you could transform your data with either projectRaster (can be slow) or gdal_warp (fast).
You need to tell the function what new projection you want the data in and it will do all the math for you. Let’s say you have sentinel data ‘r’ that are in longlat:r_proj = projectRaster(r, res = 10, crs = the proj4 string for UTM, method = ‘ngb’)What this is doing is setting the pixel size (res) as 10, and you are defining the new crs. The easiest way to get the proj4 string for UTM zone 11 is to call in a layer that already has the correct crs and copy that (you can also google which one is best for your area). It’s the same idea for vector. If you have a shapefile that isn’t in the right projection and it is an sf object:
vect_proj = st_transform (old, crs = new proj4 string)
With all of these, you must write them to file (writeraster) after you re-project. You don’t want to be projecting every time because it can take a while.
If you are trying to stack a bunch of raster layers, it is important that they have the same extent (min max xy coords). You can check this using extent(raster layer). If they are different you can actually set the extent when you project the raster. In projectRaster there is a ‘to’ argument. That allows you to project a sentinel raster for example ‘to’ the pixels of the other raster.
r_targ = raster(grid metric)
r_proj = projectRaster(sentinel raster, to = r_targ, crs = crs(r_targ))
This post was adapted from a lesson given by my supervisor.
Leave a Reply