2

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.

Picture 1enter image description here enter image description here

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)
0

1 Answer 1

3

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:

sheared

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.

3
  • 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.
    – david_p
    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.
    – Mike T
    Dec 2, 2014 at 19:56
  • But you could also use the built-in math.tan on the one coefficient.
    – Mike T
    Dec 2, 2014 at 19:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.