# How to split rasters tile size not in pixels but split in meters?

I write my raster calculator, similar in functionality to Qgis, only I will read arithmetic operations with Pyopencl I'm stay in only one place, there are two raster images that have the same size in meters, but a different pixel resolution raster1 =0.037368,-0.037368 raster2=0.0463053,-0.0463053

Since my rasters occupy more than 1 GB ,I'm forced to break them into buffers 512 to 512

``````bufferXSize = 512
bufferYSize = 512
for x in range(numX):
col = 0
print col
for y in range(numY):
xsize = rasterXSize - x*bufferXSize if x == numX-1 else bufferXSize
ysize = rasterYSize - y*bufferYSize if y == numY-1 else bufferYSize
slice = gdalData.ReadAsArray(x*bufferXSize, y*bufferYSize, xsize, ysize)
r1 = [rx, ry, rx + (geoT[1] * slice.shape[0]), ry + (geoT[5] * slice.shape[1])]
rx+=slice.shape[0]
ry+=slice.shape[1]
``````

My question is, because the two images have the same size in meters, how can I not retreat by 512 to 512 pixels, and what would the indentation of blocks be done in meters (for example, indent 25 to 25 meters on the first raster, and the same indent on the second)

• It is not so complicated. Please, see my answer. Jul 7, 2017 at 14:30

I would recommend to convert your pixel resolution into meters following this concept.

25m * 25m into pixel ratio would be 2.952.756 pixel * 2.952.756 pixel for 300dpi

300dpi resolution is 118 px/cm (Pixel per Centimeter)

Ground Resolution and Map Scale In addition to the projection, the ground resolution or map scale must be specified in order to render a map.

At the lowest level of detail (Level 1), the map is 512 x 512 pixels. At each successive level of detail, the map width and height grow by a factor of 2: Level 2 is 1024 x 1024 pixels, Level 3 is 2048 x 2048 pixels, Level 4 is 4096 x 4096 pixels, and so on.

In general, the width and height of the map (in pixels) can be calculated as:
map width = map height = 256 * 2level pixels
The ground resolution indicates the distance on the ground that’s represented by a single pixel in the map.

For example, at a ground resolution of 25 meters/pixel, each pixel represents a ground distance of 25 meters.
The ground resolution varies depending on the level of detail and the latitude at which it’s measured.

Using an earth radius of 6378137 meters (EPSG: WGS 84), the ground resolution (in meters per pixel) can be calculated as:
ground resolution = cos(latitude * pi/180) * earth circumference / map width = (cos(latitude * pi/180) * 2 * pi * 6378137 meters) / (256 * 2level pixels)

Calculating Resolution Map resolution is a function of the latitude, the zoom level, and a constant value. The constant is based on the diameter of the Earth and the equations Microsoft used to set the zoom levels. At any latitude and zoom level, you can determine the scale by using the following equation:

Map resolution = 156543.04 meters/pixel * cos(latitude) / (2 ^ zoomlevel)

For further information:

Answer is because it is not so complicated. If you read 'ReadAsArray' method documentation and first do a test with a small raster, you can get that following code works as expected.

``````from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')
filename = "/home/zeito/pyqgis_data/aleatorio1.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

print cols, rows

bufferXSize = 4
bufferYSize = 5

for i in range(0, rows, bufferYSize):
for j in range(0, cols, bufferXSize):
print "verts:", i, j
data = band.ReadAsArray(j, i, bufferXSize, bufferYSize)
print data
``````

I tried it out with a small raster (20 x 20) of next image:

After running the code at Python Console of QGIS, I got 20 matrices and their respective indexes of (row, column) where they begin.

``````20 20
verts: 0 0
[[ 1  6  4  9]
[10  8  1  2]
[ 4  5  6  9]
[ 3  8  5  3]
[ 7  9  7  7]]
verts: 0 4
[[ 4  7  3  2]
[ 4  1  1  8]
[ 5  4  1 10]
[ 2  3 10  3]
[ 6  3  7  7]]
verts: 0 8
[[ 7  9  4  8]
[ 8  1  4  5]
[ 1  4  9  8]
[ 7  2  3  8]
[ 5 10  3  1]]
verts: 0 12
[[ 5  5  3  6]
[ 7 10  8  6]
[10  6  1  9]
[ 3  7  7  7]
[ 2  2  1  3]]
verts: 0 16
[[ 4  7  9  6]
[ 3  8  8  1]
[ 5  7  9  7]
[10  6  5  3]
[10  2  5  4]]
verts: 5 0
[[7 9 1 6]
[5 4 1 8]
[7 8 8 4]
[3 9 6 9]
[1 8 9 6]]
verts: 5 4
[[ 3  4  3  9]
[ 2  4  7  4]
[ 9  7  5  3]
[10  1  9  5]
[ 2  9  8  5]]
verts: 5 8
[[ 8  4  5  9]
[10  9  3  8]
[ 4  5  7  2]
[ 1  3  7  4]
[10  5  7  3]]
verts: 5 12
[[8 1 4 1]
[8 3 2 8]
[8 6 5 2]
[9 5 2 4]
[5 4 6 7]]
verts: 5 16
[[ 5  6  6  6]
[ 1  2  4  6]
[ 4 10  6  5]
[ 5  4  4  9]
[ 1  2  8  4]]
verts: 10 0
[[ 3  6 10  3]
[10  8  6  8]
[ 7  4  9  2]
[ 4  8  9  7]
[ 3  5  9  4]]
verts: 10 4
[[ 2  9 10  6]
[ 6  4  8  4]
[ 2 10  9 10]
[ 5  4 10 10]
[ 2  4  8  8]]
verts: 10 8
[[ 7  7  8  3]
[ 3  2 10  2]
[ 7  4  3  3]
[ 1  6  9  2]
[ 7  1  2  4]]
verts: 10 12
[[ 7  7  1  5]
[10  3  1  4]
[10  4  6  6]
[ 2  7  9  6]
[ 1  4 10  7]]
verts: 10 16
[[8 1 4 9]
[9 1 2 8]
[7 3 1 7]
[4 3 8 5]
[6 3 7 3]]
verts: 15 0
[[ 2  4  8  1]
[ 7  5  2  9]
[ 7  4 10  4]
[ 9 10  1 10]
[ 3  5  4  3]]
verts: 15 4
[[10  5  2 10]
[ 1  1  2  8]
[ 2  6  1 10]
[ 9 10  4  6]
[ 1  7 10  5]]
verts: 15 8
[[ 7  9  7  8]
[10 10  6  7]
[ 7  4  9  3]
[10  2  6  7]
[ 2  1  7  3]]
verts: 15 12
[[ 3  5  5  4]
[ 8  6  5  7]
[ 4 10  1  9]
[ 8  8  9  4]
[10  6  6  1]]
verts: 15 16
[[ 4 10  2  9]
[ 4  3 10  9]
[ 3  4  9  8]
[ 5  4  9  8]
[ 1 10  6  6]]
``````

With Value Tool plugin I corroborated several matrices and all values matched as expected.

Editing Note

Complete code for producing (in this case) all 20 raster:

``````from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')
filename = "/home/zeito/pyqgis_data/aleatorio1.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

transform = dataset.GetGeoTransform()

xinit = transform[0]
yinit = transform[3]

xsize = transform[1]
ysize = transform[5]

cols = dataset.RasterXSize
rows = dataset.RasterYSize

print cols, rows

bufferXSize = 4
bufferYSize = 5

# Create gtif file
driver = gdal.GetDriverByName("GTiff")

path = "/home/zeito/pyqgis_data/"

wkt = dataset.GetProjection()

# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)

x = xinit
y = yinit

k = 1

for i in range(0, rows, bufferYSize):
for j in range(0, cols, bufferXSize):
print "verts:", i, j, k
data = band.ReadAsArray(j, i, bufferXSize, bufferYSize)
#        print data
x = xinit + j*xsize
y = yinit + i*ysize
print x, y
new_transform = (x, transform[1], transform[2], y, transform[4], transform[5])
output_file = path + 'new_raster' + str(k) + '.tif'
print output_file

dst_ds = driver.Create(output_file,
bufferXSize,
bufferYSize ,
1,
gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( data )

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(new_transform)

dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset
dst_ds = None

k += 1

dataset = None
``````

After running the code at Python Console of QGIS, I loaded only 19 first produced raster to corroborate that split was produced as expected.

• Each matrix can easily be convert in raster with numpy (to get an array) and GDAL methods. Jul 7, 2017 at 18:04