I need to derive a raster representing aspect for North America (from STRM data, originally in geographic lat/long). In most projections, true north is only map north along eg the central meridian so aspect values are distorted. I attempted calculating aspect in Mercator and then reprojectiong to eg Lambert Conformal Conic, but the aspects values are very pixelated and distorted.
The following answer by whuber suggests how slope could be calculated accurately across a continent: "Compute the slope (as rise:run or a percent) using a Mercator projection, then divide the result by the cosine of the latitude". Using SRTM Global DEM for Slope calculation?
How would this be extended to calculate aspect?
Would one use flow direction based on the corrected slope?
Or can other methods be used?
I see in Compute global DEM topographic indices and compromise projections that @whuber states "Having done this [corrected slope calculation], perform the followup calculation of flow directions, aspects, etc. using the corrected slopes."
I'm just not sure how to derive aspect from a slope raster (vs from a DEM) in ArcGIS or otherwise.
Answer
There are many possible solutions, but a convenient and accessible way will be to use an appropriate projection. What properties must it have?
It should be conformal (or very nearly so throughout the region of interest). This means there is no relative scale distortion as the bearing is varied around any fixed point but--necessarily--there will be distortion of scale from one point to the next. However, this latter distortion is going to be so small over the neighborhoods in which aspect is computed that it will have negligible influence on the accuracy. The lack of relative local distortion means any aspect computed in the projected coordinates will not be affected by the projection. (The slope, however, might be seriously wrong, which is why for geographically extensive regions it often needs to be computed separately as described elsewhere.)
It should be cylindrical. This means that true north always points in the same direction across the map. When using this common direction to designate the zero angle, aspects (which are computed in the projected coordinate system) will always be correct--or at least will have a constant offset from north, which is easily corrected afterwards.
The family of Mercator projections is conformal and cylindrical. (In fact, this more or less completely determines the Mercator projection, up to scale.) The solution for North America, then, which is covered by the Mercator, is to reproject the original DEM into Mercator and compute its aspect.
For those interested in working near poles, where the Mercator breaks down, it might be most convenient to use a Stereographic projection centered at the pole. Because this is not cylindrical, the aspect will have to be adjusted by the apparent direction of each meridian, which varies systematically from point to point in a simple way: the projected meridian always points directly towards where the pole is projected. For instance, suppose you were to use a Northern Stereographic projection with the Prime Meridian directed to the right. Then a calculated aspect of 270 degrees, which ordinarily would be interpreted as due West, actually corresponds to true North. To make this adjustment, 90 degrees has to be added to the aspect. At another location on the map, corresponding (say) to -90 degrees longitude, the meridian would point straight up to the pole: now, no adjustment is necessary. In general, in this situation, 90 degrees plus the longitude has to be added to the computed aspect (and then reduced modulo 360). But that's a simple sequence of grid operations that can be easily performed after an aspect computation.
For worldwide calculations, perform computations on overlapping areas covering the two poles and the near-Equatorial regions. Reproject those results to a common coordinate system and mosaic them. As a double-check, the aspects computed in the overlapping regions should be in close agreement throughout.
Because the projection used to compute the aspect might not be the one desired for further analysis, you might need to reproject the aspect grid. This must be done by computing a cosine and sine grid of the aspects, reprojecting them, and recombining those via an arctangent calculation. See https://gis.stackexchange.com/a/133709 for further discussion.
No comments:
Post a Comment