<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
Update<br>
<br>
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.<br>
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.<br>
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.<br>
<br>
Je note aussi qu'avec Mapnik les performances sont bien meilleures
quand on se passe du compositing et avec un tiff "TILED=YES".<br>
<br>
Alors si ça peut servir :<br>
<blockquote>#!/usr/bin/python<br>
# read a file block by block, and write band 1 into band 2, fill
band 1 with 0s.<br>
<br>
import gdal<br>
import numpy<br>
import sys<br>
<br>
from gdalconst import *<br>
<br>
ds = gdal.Open( 'DEM/3857-tiled.tif', GA_ReadOnly)<br>
<br>
band = ds.GetRasterBand(1) <br>
<br>
# Process by blocks, unless the dataset can fit into memory<br>
# Use blocks anyway, they are easy to read, and tiff is compressed
block by block<br>
<br>
block_sizes = band.GetBlockSize() <br>
x_block_size = block_sizes[0] <br>
y_block_size = block_sizes[1] <br>
<br>
xsize = band.XSize <br>
ysize = band.YSize <br>
<br>
raster = gdal.GetDriverByName('GTiff')<br>
dst_ds = raster.Create('DEM/out.tif', xsize, ysize, 2,
gdal.GDT_Byte,<br>
options=[ 'TILED=YES', 'COMPRESS=DEFLATE',
"BIGTIFF=YES", "ALPHA=YES" ])<br>
projection = ds.GetProjection()<br>
#~ print 'proj: ', projection<br>
geotransform = ds.GetGeoTransform()<br>
#~ print 'geo: ', geotransform<br>
<br>
dst_ds.SetGeoTransform( geotransform )<br>
dst_ds.SetProjection(projection)<br>
#~ dst_ds.GetRasterBand(1).Fill(0)<br>
<br>
for i in range(0, ysize, y_block_size): <br>
sys.stdout.write(str((ysize,i,100*(1-(ysize-i)/ysize),"%
left"))+'\n')<br>
sys.stdout.flush()<br>
if i + y_block_size < ysize: <br>
rows = y_block_size <br>
else: <br>
rows = ysize - i <br>
for j in range(0, xsize, x_block_size):<br>
if j + x_block_size < xsize: <br>
cols = x_block_size <br>
else: <br>
cols = xsize -j<br>
<br>
data=ds.GetRasterBand(1).ReadAsArray(j, i, cols, rows)<br>
inv=numpy.ones(data.shape, numpy.uint8) * 255<br>
data=inv-data<br>
dst_ds.GetRasterBand(2).WriteArray(data,j,i)<br>
z=numpy.zeros(data.shape, numpy.uint8)<br>
dst_ds.GetRasterBand(1).WriteArray(z,j,i)<br>
<br>
# Once we're done, close properly the dataset<br>
dst_ds = None<br>
ds = None<br>
</blockquote>
Yves<br>
</body>
</html>