Update
OK, j'ai finalement réussi à convertir un hillshade global depuis 1
bande en niveau de gris vers une bande toute noire et le hillshade dans
une bande alpha.
Sauf que, au bout de plus de 10 jours, Gdaldem y serait encore, et si
j'en crois sa barre de progression, pour quelques mois.
Par contre, avec le code ci-dessous, en 2 heures c'est plié sur un Tiff
de 100Go, donc je le passe ici si ça vous donne des idées.
Je note aussi qu'avec Mapnik les performances sont bien meilleures quand
on se passe du compositing et avec un tiff "TILED=YES".
Alors si ça peut servir :
#!/usr/bin/python
# read a file block by block, and write band 1 into band 2, fill
band 1 with 0s.
import gdal
import numpy
import sys
from gdalconst import *
ds = gdal.Open( 'DEM/3857-tiled.tif', GA_ReadOnly)
band = ds.GetRasterBand(1)
# Process by blocks, unless the dataset can fit into memory
# Use blocks anyway, they are easy to read, and tiff is compressed
block by block
block_sizes = band.GetBlockSize()
x_block_size = block_sizes[0]
y_block_size = block_sizes[1]
xsize = band.XSize
ysize = band.YSize
raster = gdal.GetDriverByName('GTiff')
dst_ds = raster.Create('DEM/out.tif', xsize, ysize, 2, gdal.GDT_Byte,
options=[ 'TILED=YES', 'COMPRESS=DEFLATE',
"BIGTIFF=YES", "ALPHA=YES" ])
projection = ds.GetProjection()
#~ print 'proj: ', projection
geotransform = ds.GetGeoTransform()
#~ print 'geo: ', geotransform
dst_ds.SetGeoTransform( geotransform )
dst_ds.SetProjection(projection)
#~ dst_ds.GetRasterBand(1).Fill(0)
for i in range(0, ysize, y_block_size):
sys.stdout.write(str((ysize,i,100*(1-(ysize-i)/ysize),"%
left"))+'\n')
sys.stdout.flush()
if i + y_block_size < ysize:
rows = y_block_size
else:
rows = ysize - i
for j in range(0, xsize, x_block_size):
if j + x_block_size < xsize:
cols = x_block_size
else:
cols = xsize -j
data=ds.GetRasterBand(1).ReadAsArray(j, i, cols, rows)
inv=numpy.ones(data.shape, numpy.uint8) * 255
data=inv-data
dst_ds.GetRasterBand(2).WriteArray(data,j,i)
z=numpy.zeros(data.shape, numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(z,j,i)
# Once we're done, close properly the dataset
dst_ds = None
ds = None
Yves
_______________________________________________
dev-fr mailing list
[email protected]
https://lists.openstreetmap.org/listinfo/dev-fr
