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.
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 ?