[OSM-dev-fr] [OSM-talk-fr] Ombrage et relief...

[yvecai]

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
-------------- section suivante --------------
Une pièce jointe HTML a été nettoyée...
URL: <http://lists.openstreetmap.org/pipermail/dev-fr/attachments/20140123/be283fed/attachment.html>