GDAL raster is rotated/flipped incorrectly

Question

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). But the raster, when opened in Qgis, looks like this:

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:

What exactly am I doing wrong ?

Answer 1

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: