[gdal-dev] how to create just the msk file from a rgba vrt file

[Duarte Carreira]

Hi Even, and many thanks! Your pointers have solved my question.

As for the problem building overviews it seems that it occurs when you enter a <OVERVIEW> tag pointing to a "regular" image, that doesn't have the INTERNAL_MASK_FLAGS_# set.

On the other hand these flags allowed me to create a mask file even faster than using gdal_translate.

The process is:
1) get a shapefile that represents your mask (maybe use gdaltindex and dissolve the result)
2) use gdalinfo to get extent and cell size of your vrt mosaic 
3) use gdal_rasterize to get your mask raster, burning the shapefile to a new image initiated with 0 and burnt with 255
4) set the metadata tags to make this raster a mask (I used a small python script for this)
5) rename your raster mask to match your vrt mosaic (eg mosaic.vrt.msk)

That's it. Instant (almost) mask file: instead of 4h41min it took 10min to create my mask. (The vrt mosaic references 487 rgb tiff files, totaling 240001x250000 pixels in size.)

The command to build the mask was the following:

gdal_rasterize --config GDAL_CACHEMAX 512 -of GTiff -co compress=deflate -co tiled=yes -l quads_diss_etrs -te -79997.180 -229999.390 40002.909 -104999.818 -tr 0.499994258518990 0.499994258518990 -ot Byte -burn 255 -init 0 -co INTERLEAVE=BAND -co NBITS=1 quads_diss_etrs.shp maskfile.tif

I'm not sure about the NBITS command effect though...

Then to insert the mask flags I used this python script:

#################################################
import gdal
from gdalconst import *
import sys

gdal.UseExceptions()
try:
  filename = sys.argv[1]

  dataset = gdal.Open(filename, GA_Update)
  dataset.SetMetadataItem( "INTERNAL_MASK_FLAGS_1", "2", None )
  dataset.SetMetadataItem( "INTERNAL_MASK_FLAGS_2", "2", None )
  dataset.SetMetadataItem( "INTERNAL_MASK_FLAGS_3", "2", None )

except Exception, e:
  print "Error: " + str(e)
  
  
dataset = None
##################################################

Thanks again.
Duarte



-----Mensagem original-----
De: Even Rouault [mailto:even.rouault at mines-paris.org] 
Enviada: domingo, 9 de Fevereiro de 2014 09:56
Para: gdal-dev at lists.osgeo.org
Cc: Duarte Carreira; Brian Case
Assunto: Re: [gdal-dev] how to create just the msk file from a rgba vrt file

Le vendredi 07 février 2014 10:55:06, Duarte Carreira a écrit :
> Thanks Brian.
> 
> But this way you rewrite the whole image to disk. It uses lots of disk 
> space and takes forever.
> 
> I want to avoid that and just get the mask out to a msk file as fast 
> as possible.
> 
> I don't want to convert my rgba vrt mosaic.
> 
> The final objective is to get a msk file I can use with the simple rgb 
> vrt mosaic. Then I'll be able to build overviews with jpeg/ycbcr 
> compression which I can't do with rgba vrt because of the 4 bands.
> 
> For now I have tried 3 ways:
> 
> 1) use gdal_rasterize to create a mask directly from the mask polygon 
> shapefile. Then just edit the rgb vrt mosaic and add a  maskband to it.
> The problem here is gdaladdo does not honor the maskband. This is the 
> fastest way I know, pitty it doesn't work in the end.

I'd be curious that you provide ways of reproducing this. The overview computation code has explicit code to deal with mask bands.

> 
> 2) use gdal_translate like you suggested but use -scale to write all 
> 0s in all 3 bands, and compress with deflate. You get a valid mask and 
> a very, very small useless mosaic. This works but takes a while still.
> 
> 3) use gdal_translate like you suggested but exaggerate the jpeg 
> compression so it errors out (jpeg_quality=15). You get an invalid 1kb 
> mosaic and an apparently good msk. But it's corrupted in some way. So 
> doesn't work.
> 
> I think #1 has potential. If there was a way to somehow turn the tif 
> created by gdal_rasterize into a "true" mask file and have it honored 
> by gdaladdo we would have a winner.
> 
> Maybe there's a way to directly export the alpha band from the rgba 
> vrt mosaic to a mks file without writing anything else? That I guess 
> would be the fastest way of all.

You can generate a valid .msk file with the following gdal_translate command :

gdal_translate -b 4 rgba.tif out.tif.msk \ -mo "INTERNAL_MASK_FLAGS_1=2" -mo "INTERNAL_MASK_FLAGS_2=2" \ -mo "INTERNAL_MASK_FLAGS_3=2" -co COMPRESS=DEFLATE -co INTERLEAVE=BAND \ -co NBITS=1

You need to provide as many -mo "INTERNAL_MASK_FLAGS_X=2" option as there are bands (so 3 for a RGB dataset as in the above example). This is the important option that will make a .msk file being recognized as a mask band.

Even

--
Geospatial professional services
http://even.rouault.free.fr/services.html