1

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

enter image description here

enter image description here

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. – xunilk Jul 7 '17 at 14:30
1

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)

(Source: https://www.blitzrechner.de/pixel-zentimeter-umrechnen/)

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:

https://msdn.microsoft.com/en-us/library/bb259689.aspx
and
https://msdn.microsoft.com/en-us/library/aa940990.aspx

0

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:

enter image description here

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.

enter image description here

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.

enter image description here

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

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.