# Slanted raster along the line? Wind erosion model and windbreak effectivity

I'm working on wind erosion modelling and windbreak potential. Imagine windbreak, which has some values of optical porosity (in %) and this values determinates the effectivity of windbreak in specific part of windbreak. I want to compute this effectivity in some distance behind this windbreak and make a raster from it. That's not a problem, if the windbreak is oriented orthogonal and wind direction is upright on the windbreak (picture 1).

The main problem is when windbreak (line) do not lays oriented orthogonal and / or wind direction isn't upright on the windbreak (picture 2 and 3). I can't find any desription for module GDAL about this, specialy for SetGeoTransform function, where is 3rd and 5th argument define rotation of raster cells. My idea is that y axis cell rotation is usefull when windbreak (line) is not orhtogonal oriented, so whit this I can snap raster along the line, and x axis cell rotation can be used for turn a raster by the wind direction.

I hope that you can understand what I mean. Anyway I can explain more about this when it's necessary.

Here's my code what i have done.

``````import numpy
import os
import arcpy
from osgeo import gdal, gdal_array, osr

os.chdir("d:\wind")

# matrix contains values of optic porosity for windbreak, 6 horisontal layers, each layer 9 vertical values
porosity_matrix = numpy.matrix(
[[96, 95, 96, 78, 95, 100, 93, 70, 74], [73, 57, 45, 42, 71, 80, 65, 43, 50],
[47, 29, 32, 30, 31, 30, 29, 34, 35], [41, 48, 32, 40, 61, 52, 37, 53, 47],
[28, 31, 28, 29, 57, 49, 36, 46, 34], [28, 48, 37, 24, 39, 27, 30, 24, 20]])

# definition of line coordinates, finds a first point of raster extent
arcpy.env.workspace = r"d:\wind"
layer = "windbreak.shp"
desc_layer = arcpy.Describe(layer).shapeFieldName
cur_layer = arcpy.SearchCursor(layer)
for x in cur_layer:
obj = x.getValue(desc_layer)
FP = obj.firstPoint
LP = obj.lastPoint
coordFP = [FP.X, FP.Y] # coordinates of line's first point
coordLP = [LP.X, LP.Y] # coordinates of line's last point
coordDiference = [(coordFP[0] - coordLP[0]), (coordFP[1] - coordrLP[1])]
y_res = int(coordDiference[1]) / numpy.shape(prosity_matrix)[1] # definition of raster's cells y resolution (based on measured data)

# function for compute effectivity and make raster layer
def effectivity(in_matrix, row_number, ras_name):
for i in range(9):
Y = in_matrix[row_number, i]
values_row = []
for X in range(5, 405, 5): # X is a distance, 5m step
Z = 46.1894 + 0.1709 * X - 0.0008 * X * X - 0.0004 * X * Y - 0.000094226 * Y * X # windbreak effectivity empiric equation
Z100 = 100 - Z # Z is a effectivity of windbreak for one value of porosity
if Z100 > 100:
Z100 = 100
values_row.append(Z100)
if i == 0:
out_matrix = numpy.matrix((values_row))
else:
matrix_temp = numpy.matrix((values_row))
out_matrix = numpy.concatenate((out_matrix, matrix_temp))
## create raster
name = str(ras_name) + ".tif"
nrows, ncols = numpy.shape(out_matrix)
xres, yres, xmin, ymax = (5, y_res, FP.X, FP.Y)
# important part
geotransform = (xmin, xres, 0, ymax, 0, -yres) # this is important for inclinate raster
ras_name = gdal.GetDriverByName("GTiff").Create(name, ncols, nrows, 1, gdal.GDT_Float64)
ras_name.SetGeoTransform(geotransform)
# important part ends
srs = osr.SpatialReference()
srs.ImportFromEPSG(5514)
ras_name.SetProjection(srs.ExportToWkt())
ras_name.GetRasterBand(1).WriteArray(out_matrix)
ras_name = None

# for every horisontal layer of optical porosity creates a raster
for q in range(6):
name = "windbreak_" + str(q)
effectivity(prosity_matrix, q, name)
print "Created raster of effectivity no.: " + str(q)
``````

You need to apply a shear to your affine transformation matrix used for SetGeoTransform. One package that is intended to work with these is `affine`.

For example:

``````from affine import Affine

xmin = 1; xres = 2; ymax = 3; yres = 4  # dummy values

# construct geotransform matrix
geotransform = Affine(
xres, 0, xmin,
0, -yres, ymax)

# or import from raster using: geotransform = Affine.from_gdal(*ras_name.GetGeoTransform())

print(repr(geotransform))
# Affine(2.0, 0.0, 1.0,
#        0.0, -4.0, 3.0)

# apply shear, e.g. 10 degrees
geotransform *= Affine.shear(10, 0)

# also works: geotransform = geotransform * Affine.shear(10, 0)

print(repr(geotransform))
# Affine(2.0, 0.0, 1.0,
#        -0.7053079228338599, -4.0, 3.0)

# reorganise coefficients for GDAL
ras_name.SetGeoTransform(geotransform.to_gdal())
``````

The two rasters (regular and sheared) are shown here, where the yellow square is the sheared result:

In summary, you don't need to use additional packages like `affine`, but you do need to understand what's happening with the affine transform matrices.

• Thanks! That exactly what I need. And I assume that the same way I could inclinate raster along the line. So I have to make transformation two times. First for "snap" raster along the line and second for "shear" raster by the wind direction. But can I ask you a dumb question? Where can I get the affine modul for python? I found this, but in the file there wasn't setup.py, so I don't know how to install it into python. I'll try to search more. Dec 2, 2014 at 14:46
• For Windows, I've installed pip, then open a command prompt, `cd c:\Python27\Scripts` (or wherever you have Python installed), then `pip install affine`. Alternatively you could extract the zip file, then CD into the folder, then run `python setup.py install`. Dec 2, 2014 at 19:56
• But you could also use the built-in `math.tan` on the one coefficient. Dec 2, 2014 at 19:59