In [80]:

    

import osgeo.gdal
import osgeo.osr

try:
    from PIL import Image
except ImportError:
    import Image


    

In [2]:

    

drivers = {osgeo.gdal.GetDriver(i).ShortName:osgeo.gdal.GetDriver(i) for i in range(osgeo.gdal.GetDriverCount())}


    

In [3]:

    

memdriver = drivers['MEM']
tiffdriver = drivers['GTiff']


    

In [20]:

    

src_ds = memdriver.Create("", 1000, 1000, 1, eType=osgeo.gdal.GDT_Int32)
src_ds.RasterCount
geotransform = (255580.19380886,
 4.7338744607343,
 0.0,
 6247234.441056012,
 0.0,
 -4.7338744607343)
src_ds.SetGeoTransform(geotransform)
src_srs = osgeo.osr.SpatialReference()
epsg = 22234
src_srs.ImportFromEPSG(epsg)
src_ds.SetProjection(src_srs.ExportToWkt())


    

    
Out[20]:





0

In [20]:

    




    

In [21]:

    

band = src_ds.GetRasterBand(1)
m = band.XSize
n = band.YSize
# fill with indices
arr = np.arange(n*m).reshape(m,n)
band.WriteArray(arr)
arr
plt.imshow(arr, vmin=arr.min(), vmax=arr.max())


    

    
Out[21]:





<matplotlib.image.AxesImage at 0x12951490>






    

In [61]:

    

srs = 900913
width=2573
height=1706

grayscale_src = (src_ds.RasterCount == 1)

# Prepare output gdal datasource -----------------------------------

%timeit area_ds = tiffdriver.Create('/vsimem/output.tiff', width, height, 1, eType = osgeo.gdal.GDT_Int32)

if area_ds is None:
    raise Exception('uh oh.')


merc = osgeo.osr.SpatialReference()
merc.ImportFromEPSG(epsg)

area_ds.SetProjection(merc.ExportToWkt())


gtx = [x, w/width, 0, y, 0, h/height]
gtk = list(geotransform)
gtx[0] += 3
gtx[1] += 0.1
gtx[3] += 0.3
gtx[5] -= 0.21
area_ds.SetGeoTransform(gtx)
band = area_ds.GetRasterBand(1)
band.SetNoDataValue(-9999)
band.WriteArray(np.zeros((band.YSize, band.XSize))-9999)


# Adjust resampling method -----------------------------------------

resample = osgeo.gdal.GRA_NearestNeighbour

# Create rendered area ---------------------------------------------

%timeit osgeo.gdal.ReprojectImage(src_ds, area_ds, src_ds.GetProjection(), area_ds.GetProjection(), resample, 1000000, 0.1 )
area = area_ds.GetRasterBand(1).ReadAsArray()
area = np.ma.masked_equal(area, -9999)
tiffdriver.Delete('/vsimem/output.tiff')


    

    




10 loops, best of 3: 33.3 ms per loop
100 loops, best of 3: 9.14 ms per loop






    
Out[61]:





0

In [62]:

    

%%timeit
B = np.ma.take(arr, area)


    

    




1 loops, best of 3: 202 ms per loop

In [79]:

    

?np.take


    

In [ ]: