2

I created a GeoTIFF file with size 3600x3600, I verified that the file I have on harddrive is indeed 3600x3600. But the raster size that QGIS displays is 4554x2277.

This makes the raster sum wrong, (407Million instead of the correct 509Million). Why does QGIS show a different raster size here? Why the particular number 4554?

Script I used for creating the raster file.

from osgeo import gdal
from osgeo import osr
import numpy
import math

fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)

xdim = 3600
ydim = 3600

dst_filename = "test.tif"
dst_ds = driver.Create(dst_filename, xsize=xdim, ysize=ydim,bands=1, eType=gdal.GDT_Float32)


x0 = -180
y0 = -90
x_per_pixel = 360.0 / xdim
y_per_pixel = 180.0 / ydim

dst_ds.SetGeoTransform([x0, x_per_pixel, 0, y0, 0, y_per_pixel])

srs = osr.SpatialReference()
#srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("WGS84")
proj_str = srs.ExportToWkt()
print('proj: ',proj_str)
dst_ds.SetProjection(proj_str)

raster = numpy.zeros((xdim, ydim), dtype=numpy.float)
earth_a =  6378137 / 1000
earth_b = 6356752.3142 / 1000
for i in range(xdim):
    for j in range(ydim):
        lat = y0 + y_per_pixel * (j + 0.5)
        a = 2 * math.pi * earth_a / xdim
        b = math.pi * earth_b / ydim
        raster[j][i] = math.cos(math.radians(lat)) * a * b


dst_ds.GetRasterBand(1).WriteArray(raster)
# Once we're done, close properly the dataset
dst_ds = None

What QGIS shows in Properties->Information of the raster layer

Name    test
Path    I:\gis\worlddiv\py\test.tif
CRS EPSG:4326 - WGS 84 - Geographic
Extent  -180.0000000000000000,-90.0126558050849894 : 180.0253116101699788,90.0000000000000000
Unit    degrees
Width   4554
Height  2277
Data type   Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF
Dataset Description 
Compression 
Band 1  
STATISTICS_APPROXIMATE=YES
STATISTICS_MAXIMUM=61.75239944458
STATISTICS_MEAN=38.018121292831
STATISTICS_MINIMUM=0.026944568380713
STATISTICS_STDDEV=19.920540978552
STATISTICS_VALID_PERCENT=100
More information    
Dimensions  X: 4554 Y: 2277 Bands: 1
Origin  -180,90
Pixel Size  0.07905694150420948529,-0.07905694150420948529
2
  • The two answers provided so far are great. Please allow me to add another possibility. I'm doing it as a comment rather than answer because I'm not confident that my idea is correct, but in the spirit of helpfulness, here goes: I note that your Properties > Information shows "STATISTICS_APPROXIMATE=YES". Perhaps this approximation is the source of your difference. I previously asked about such approximations, and received some great info. Maybe it will assist you: gis.stackexchange.com/questions/321408/…
    – Stu Smith
    Dec 20, 2020 at 22:11
  • @StuSmith Nice idea! Butt the spatial dimensions of the raster file are not part of those statistics, they are used for the actual raster values. Dec 20, 2020 at 22:46

2 Answers 2

2

It seems like QGIS is warping the original image, that is created as South-up, into a North-up image. The code makes the same thing than running gdalwarp with default settings. Gdalwarp is not guaranteed to keep the resolution and by default it creates square pixels with equal x and y size. Therefore the warped image will be two times wider that high also as pixels, not only as geographic units (360 vs. 180 degrees). The x_per_pixel and y_per_pixel settings which were used for creating the original image gets overrun.

A test with gdalwarp

gdalwarp test.tif testwarped.tif
Creating output file that is 4554P x 2277L.

Edit your script a bit and place origin at the top left corner and make the y in the geotransform negative. That way a North-up image is created and QGIS will accept it without warping.

x0 = -180
y0 = 90
x_per_pixel = 360.0 / xdim
y_per_pixel = -180.0 / ydim

This way you'll create a common north-up image and QGIS seems to get it right.

enter image description here

The QGIS behavior with your original image with origin at lower left corner is indeed odd and properties show different values than than those reported by gdalinfo. This may be a bug in QGIS. Probably your original image is valid even it is of uncommon sort but I am not sure about that.

gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 3600, 3600
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000000,-90.000000000000000)
Pixel Size = (0.100000000000000,0.050000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Lower Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Upper Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Lower Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=3600x1 Type=Float32, ColorInterp=Gray
2

I think it may be that when creating the raster more pixels are generated, perhaps due to some inclination due to georeferencing. Load the raster into QGIS and check its properties in the python console using both PyQGIS and Gdal.

Run this code in the console, on the active raster layer, it compares the dimensions with PyQGIS and GDAL, they must be equal.

At the end it shows you the rotations by the X and Y axes, if they are different from zero (0), it is the cause of a greater number of pixels, that must be non data

import gdal

layer=iface.activeLayer()
#With PyQGIS
#rows
r=layer.height()
#columns
c=layer.width()
#x resolution (pixel)
px=layer.rasterUnitsPerPixelX()
#y resolution (pixel)
py=layer.rasterUnitsPerPixelY()
print('Dimensions PyQGIS')
print('rows: ',r,', cols: ',c,', resolutionX: ',px,', resolutionY: ',py)

#With Gdal
glayer=gdal.Open(layer.source())
#columns
gc= glayer.RasterXSize
#rows
gr= glayer.RasterYSize

geot= glayer.GetGeoTransform()
#x resolution (pixel)
gpx=geot[1]
#y resolution (pixel)
gpy=geot[5]
print('Dimensions Gdal')
print('rows: ',gr,', cols: ',gc,', resolutionX: ',gpx,', resolutionY: ',gpy)
print('Rotations: ',geot[2],' , ',geot[4])

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.