3

I am producing a raster file in python. However, the rasters are incorrect: they appear to be rotated 90° clockwise. I have provided sample code and results to illustrate the problem.

This is basically how I wanted my raster to look like (using his2d). Desired Output But the raster, when opened in Qgis, looks like this: rotated raster

As you can see, it seems to be rotated 90 degrees clockwise. Below is an example of the code I am using. Note I have just created random data because the actual data is much larger. I am not using a specified coordinate system.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from osgeo import gdal

# X-coordinates (in meters)
dataX = np.random.uniform(-50,50,100)
# Y-coordinates (in meters)
dataY = np.random.uniform(-50,50,100)
# some other variable
dataZ = np.random.uniform(-50,50,100)

# cell dimensions
binWidth = 4
binLength = 20

# min/max coordinates
xMin = min(dataX)
xMax = max(dataX)
yMin = min(dataY)
yMax = max(dataY)

# create histogram to compare with raster
plt.hist2d(dataX, dataY, bins = [np.arange(xMin, xMax+binWidth, binWidth), np.arange(yMin, yMax+binLength, binLength)], weights = dataZ)

# binning
statistic,x_edge,y_edge,binnumber = stats.binned_statistic_2d(dataX, dataY, dataZ, 'mean',bins = [np.arange(xMin, xMax+binWidth, binWidth), np.arange(yMin, yMax+binLength, binLength)] )

# get number of rows and columns
nrows,ncols = np.shape(statistic)

# create the raster  
geotransform=(xMin,20,0,yMax,0, -4)  
output_raster = gdal.GetDriverByName('GTiff').Create('testRaster.tif',ncols, nrows, 1 ,gdal.GDT_Float32) 
output_raster.SetGeoTransform(geotransform)  
output_raster.GetRasterBand(1).WriteArray(statistic)  
output_raster.FlushCache()

plt.show()

I have played around with the 'SetGeoTransform' parameters (I am assuming this is where the problem is) multiple times but none of my attempts have gotten the desired results. For example, switching the pixel sizes by setting geotransform=(xMin,4,0,yMax,0, -20), yields the following result:

enter image description here

What exactly am I doing wrong ?

  • By comparing with gdal.org/gdal_tutorial.html it seems that you define the pixel size to be 20 wide and 4 high. – user30184 Jun 18 '18 at 10:26
  • I did try switching it before : geotransform=(xMin,4,0,yMax,0, -20) , but then the raster appears stretched. I have added an additional image to show it looks if I swap the pixel sizes as you suggested. – user3451660 Jun 18 '18 at 11:00
  • Please clarify a) what is the size of the raster that you want to create as count of pixels horizontally (width) and vertically (height) and b) what is the pixel size that you want to have in georeferenced units. What geotransform does is to set the latter but perhaps your aim is not to stretch the pixels. – user30184 Jun 18 '18 at 11:22
  • a) I want size of the raster dependent on whatever the X and Y values are, so in terms of pixel count , this varies: nrows,ncols = np.shape(statistic) b) I want the pixels to have a length of 20 m (Y-distance)and a width of 4 m(X-distance), as shown in the histogram. – user3451660 Jun 18 '18 at 12:00
  • Bad result was produced by a wrong reshape and it also needs a new geotransform: geotransform=(xMin,4,0,yMax,0, -20). Please, see my answer. – xunilk Jun 18 '18 at 15:26
2

Bad result was produced by a wrong reshape and it also needs a new geotransform: geotransform=(xMin,4,0,yMax,0, -20). Following code (with my own path) works as expected:

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from osgeo import gdal

# X-coordinates (in meters)
dataX = np.random.uniform(-50,50,100)
# Y-coordinates (in meters)
dataY = np.random.uniform(-50,50,100)
# some other variable
dataZ = np.random.uniform(-50,50,100)

# cell dimensions
binWidth = 4
binLength = 20

# min/max coordinates
xMin = min(dataX)
xMax = max(dataX)
yMin = min(dataY)
yMax = max(dataY)

# create histogram to compare with raster
plt.hist2d(dataX, dataY, bins = [np.arange(xMin, xMax+binWidth, binWidth), np.arange(yMin, yMax+binLength, binLength)], weights = dataZ)

# binning
statistic,x_edge,y_edge,binnumber = stats.binned_statistic_2d(dataX, dataY, dataZ, 'mean',bins = [np.arange(xMin, xMax+binWidth, binWidth), np.arange(yMin, yMax+binLength, binLength)] )

# get number of rows and columns
nrows,ncols = np.shape(statistic)

new_statistic = np.reshape(statistic, (ncols,nrows))

nrows,ncols = np.shape(new_statistic)

# create the raster  
geotransform=(xMin, 4, 0, yMax, 0, -20)  
output_raster = gdal.GetDriverByName('GTiff').Create('/home/zeito/pyqgis_data/testRaster.tif', ncols, nrows, 1 ,gdal.GDT_Float32) 
output_raster.SetGeoTransform(geotransform)  
output_raster.GetRasterBand(1).WriteArray(new_statistic)  
output_raster.FlushCache()

plt.show()

After running it at Python Console of QGIS I got desired result:

enter image description here

  • I might be being stupid, but the new raster doesn't appear to be a match at all. I've tried it with different data sets and the although the dimensions and cell size is correct, they look completely different to the histogram.Before it was an exact match to the histogram but incorrectly shaped. :( – user3451660 Jun 19 '18 at 7:29
  • I will review it with Value Tool plugin later to corroborate if data matrix and raster match. – xunilk Jun 19 '18 at 8:48
  • I think I managed to solve it ! With new_statistic = np.reshape(statistic, (ncols,nrows)), the values will be arranged at the incorrect positions. Looking at the arrays, we actually want the transpose of it, so I tried new_statistic = statistic.transpose(), but noticed that the raster was flipped about the X axis(?). So I then added new_statistic = np.flip(new_statistic,0) to flip it around and it finally matched the histogram :) – user3451660 Jun 19 '18 at 9:06
  • They never could match because original histogram has a resolution of 25x5 and expected raster has a resolution of 5x25. By using flip numpy method you fixed original histogram to match with expected raster. By the way, GDAL handles rows and columns in a different way to linear algebra convention (rows, columns). In GDAL, X values are referred as columns and Y values as rows. This is a common source of error. – xunilk Jun 19 '18 at 9:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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