lundi 1 décembre 2014

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.


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)




Aucun commentaire:

Enregistrer un commentaire