In [1]:

    

import osgeo.osr


    

In [2]:

    

import numpy as np


    

In [3]:

    

srs="EPSG:3857"
bbox=(415817.4338713588,6746026.368336517,454953.1923533691,6770486.217387771)
height=1280
width=2048


    

In [64]:

    

rc_srs = osgeo.osr.SpatialReference()
src_srs.ImportFromEPSGA(int(srs.split(':')[1]))
dst_srs = osgeo.osr.SpatialReference()
dst_srs.ImportFromEPSGA(28992)
src2dst = osgeo.osr.CoordinateTransformation(src_srs, dst_srs)


    

    




1000 loops, best of 3: 317 µs per loop

In [65]:

    

(xmin, ymin, xmax, ymax) = bbox
xmin_dst, ymin_dst, _ = src2dst.TransformPoint(xmin,ymin)
xmax_dst, ymax_dst, _ = src2dst.TransformPoint(xmax, ymax)


    

    




100000 loops, best of 3: 6.03 µs per loop

In [66]:

    

xmin_src, ymin_src = (33940.000, 384010.000) #(  3d38'47.06"E, 51d25'59.13"N)
xmax_src, ymax_src = (87955.000, 431550.000) #(  4d24'50.49"E, 51d52' 9.89"N)
dx_src, dy_src = (5, 5) #Pixel Size = (5.000000000000000,-5.000000000000000)
x_src = np.arange(xmin_src, xmax_src, dx_src)
y_src = np.arange(ymin_src, ymax_src, dy_src)

x_start, x_end = bisect.bisect(x_src, xmin_dst), bisect.bisect(x_src, xmax_dst)
y_start, y_end = bisect.bisect(y_src, ymin_dst), bisect.bisect(y_src, ymax_dst)
x_step = (x_end - x_start) // width
y_step = (y_end - y_start) // height


    

    




10000 loops, best of 3: 78.8 µs per loop

In [61]:

    

S = np.s_[x_start:x_end:x_step, y_start:y_end:y_step]