# Calculating field length in certain directions (raster)?

Currently I'm working on winderosion. For my purpose I like to calculate the length of agricultural fields in certain wind direction. I would like to do it for raster.

To be more precise, I have agricultural areas as SHP and can convert them to raster. After it I'll calculate the distance of one boarder of the area to another in a certain direction (for instance 270° West). At the end I'll weight the different rasters (with field length in different directions) to get a mean feld length.

So my idea was to designate those cells within the area "1", those outside "0" and use the raster calculator to start calculating with the first "1" and add 20 (m = the resolution of one cell). The next cell should take the value of the previous one plus 20 and so on until it reaches "0" again.

This is just the idea but I have no idea how to do it in the raster calculator of ArcGIS/QGIS/SAGA and use a certain direction for calculating.

Do you have any proposals for me?

I also thought about neighborhood or zonal statistics.

• Interesting! It is possible to calculate neighborhood statistics in ArcMap in only one direction (e.g. "all cells to the west") -- use an irregular or weighted neighborhood when running Focal Statistics tool (using SUM as the statistic of interest). The challenge would be finding a "window" size that is large enough to capture the whole field, but small enough to not pick up part of a nearby field. – Erica Jun 13 '14 at 12:37
• Thanks for your answer, Erica. I already tried Focal Statistics with an irregular kernal file but as you said the challenge is to define a "window" for calcuation. I tried just a simple one like 1 0 and thought about rotating the grid before and rotate it back afterwards. It didn't worked out. – xmisx Jun 13 '14 at 12:42
• A somewhat similar problem is discussed in this ArcGIS forum post, but is more about measuring polygons that measuring a contiguous raster area. It might work depending on the size/shape/irregularity of the polygons. – Erica Jun 13 '14 at 12:52

## 1 Answer

The answer depends on what "field length" means.

Suppose the length is the longest segment, oriented with the wind, that is contained within the field. First rotate the data so that the wind direction now appears horizontal (west to east), then apply any of the solutions to this problem for the special case posted at How to calculate the maximum distance within a polygon in x-direction (east-west direction). (There are solutions there both for the original shapefile format and solutions for the raster representation.)

Suppose the length is the average distance across the field covered by random air particles blown by the wind. Because the area of the field is the integral of these distances, the length along the wind is the area divided by the width of the field (perpendicular to the wind's direction). Areas are easy to come by with a GIS; we need to compute the width. Let `a` denote the orientation of the wind. A unit vector perpendicular to that direction is

``````v = (sin(a), -cos(a))
``````

The distance (relative to some arbitrary origin) perpendicular to the wind at a location with coordinates (x,y) is the dot product of (x,y) with v, given by

``````d(x,y) = x * sin(a) - y * cos(a)
``````

Therefore the width perpendicular to the direction equals the maximum value of d(x,y) attained the field minus the minimum value of d(x,y) in the field. A map algebra workflow to implement this solution begins with the indicator grid of the fields described in the question: call this [Fields].

1. RegionGroup [Fields] to produce a grid with a unique identifier for each contiguous field. Call the result [Zones].

2. Obtain two grids, [X] and [Y], giving projected coordinates in the desired units of measurement. (How to do this varies with the software and its version; several methods are described in other threads on this site.)

3. Compute the local linear combination `[D] = sin(a) * [X] - cos(a) * [Y]`. In this expression [X] and [Y] are grids while `sin(a)` and `cos(a)` are numbers computed from the given wind direction `a`. (Alternatively, if the wind direction varies from field to field, `a` would be a grid giving the wind orientation at every point.)

4. Find the Zonal Max and Zonal Min of [D] relative to [Zones]. Call these [DMax] and [DMin], respectively. They can be represented either as grids or, more efficiently, as attributes in a table indexed by the zones.

5. Compute the field widths `[DMax] - [DMin] + w`; call the result [Width]. In this formula `w` equals the width of one cell perpendicular to the wind. A reasonable estimate is the cellsize `c` itself. (The reason this adjustment is needed becomes clear when you contemplate a field that is, say, perfectly horizontal in the grid and just one cell high. In such cases [DMax] and [DMin] will coincide, suggesting the field's width is zero, whereas in fact it is one.) This value estimates the length of the black dashed vector in the figure: there will be one such estimate for each distinct field of the [Zones] grid.

6. Notice that [Zones] has a "Count" field giving the number of cells per field. Because the area of a field equals c^2 times the number of cells, the computation is completed by means of the ratio

``````c^2 * [Zones].Count / [Width]
``````

This will be a raster calculation when the zonal statistics are stored in grid format or a field calculation otherwise: just join the [Zones] attribute table to the zonal statistics table from step (4) and proceed.

The actual syntax of these calculations depends on what software is used (and even on what version it is), on the interface to it (e.g., on whether separate tools are invoked or the process is scripted in a language like Python or with a command line interface in ArcGIS or GRASS) and on whether raster calculations or field calculations are employed in steps (4) through (6). Consult the appropriate help pages for the details.

• Thanks for your brilliant help, whuber. It worked out perfectly and bring exactly the results I like to have! For my purpose the second approach was more suitable. I'm not very familiar with mathematics so I would like to ask you, if you could explain the equation for the the unit vector v=(sin(a), -cos(a)) a bit more detailed for me. Since I was searching hard for step 2, I'll link to the thread where you already gave a answer to that topic: forums.arcgis.com/threads/… – xmisx Jun 15 '14 at 10:37
• People define the sine and cosine of a bearing a as the components of the unit vector in the direction of a. Also, by definition, two vectors (x,y) and (u,v) are perpendicular if and only if their dot product xu + yv equals zero. It is simple to check that the dot product of (sin(a), -cos(a)) with the direction vector (cos(a), sin(a)) is zero. That's all you need to know theoretically. As a practical matter, just make sure that a is measured in the units (usually either degrees or radians) expected of your system's implementation of the sine and cosine functions. – whuber Jun 15 '14 at 19:16
• Thanks again. Do you have any idea how to solve the problem according to the idea I posted in my first question: to designate the distance for every cell to the border of the field (or to a hedge etc.) in the direction of the wind? Your solutions were more zonal and median...still very acceptable but I can't stop thinking about counting the distance for every single cell ;) – xmisx Jun 16 '14 at 19:54
• To find those distances, first rotate the field shapes so that the wind appears to blow perfectly horizontally (perfectly vertically would do, too). Set up a grid of constant flow direction values within the fields for a flow against the wind's direction, create another grid of unit rainfall rates within the fields, and compute their flow accumulation grid (multiplied by the cellsize). Unrotate the result. – whuber Jun 16 '14 at 19:58
• Increasing resolution is the best solution. For an approximate solution, buffer the polygons internally by about one-half a cellsize and rasterize that: this procedure will guarantee all polygons are separated. Afterwards, add the buffer radius to the distances-to-the-edge that you calculate. A painful but accurate solution is to write a script that loops over each of the original shapes: this will have a lot of computational overhead and (obviously) requires more scripting. – whuber Jun 18 '14 at 13:31