[GRASS-user] Landsat image tiffs

[Dylan Beaudette]

On Friday 18 August 2006 15:06, Zenon Panoussis wrote:
> This falls in the cracks between grass, gdal and mapserver. I'm
> posting it here on the assumption that this is the place where
> someone most likely might have run into the same problem and
> solved it.
>
> I did the following ($view is a landsat set, in this case 007068):
>
>  # Import RGB & pan
>  for band in 10 20 30 80
>  do
>    file=`ls -1 $satbase/$view |grep nn$band`
>    r.in.gdal input=$satbase/$view/$file output=$view\_$band
>  done
>
>  # Sharpen
>  i.fusion.brovey -l ms1=$view\_10 ms2=$view\_20 ms3=$view\_30 \
>      pan=$view\_80 outputprefix=$view
>
>  # Beautify
>  i.landsat.rgb red=$view.red green=$view.green blue=$view.blue
>
>  # Join
>  r.composite red=$view.red green=$view.green blue=$view.blue \
>      output=$view_15me
>
>  # Check
>  d.mon start=x0
>  d.rast $view\_15me
>  # yes, it looks good
>
>  # Export
>  r.out.gdal input=$view\_15me format=GTiff output=$view\_15meB.tif \
>      type=Byte createopt="TFW=YES,TILED=YES"
>  r.out.gdal input=$view\_15me format=GTiff output=$view\_15meUI.tif \
>      type=UInt16 createopt="TFW=YES,TILED=YES"
>  r.out.tiff -t -l -v input=$view\_15me output=$view\_15meROT.tif
>
> Both r.out.gdal commands cause
> http://grass.itc.it/pipermail/grassuser/2006-August/035750.html
> but continue and produce the tiffs. r.out.tiff does not issue the
> same warning.
>
> My mapfile says
>
>  LAYER
>    NAME "landsat"
>    TYPE RASTER
>    STATUS DEFAULT
>    DATA "007068_15me[suffix].tif"
>    OFFSITE 0 0 0
>  END
>
> Using any one of the three tiffs, I get
>
>   shp2img -m world.map -o /dev/null -layer_debug landsat 10
>   [Fri Aug 18 16:21:39 2006].381669 msDrawRasterLayerLow(landsat):
> entering.
>
> That's it. mapserver is happy, debug shows no errors, there's only
> one little problem: the image doesn't show on the map. It's simply
> not there. mapserver reads it fine, but looks at the wrong band.
> Or it doesn't know where to find the colour table. Or ignores the
> tiff's colour table and, since it doesn't have one of its own, throws
> away the data. Removing OFFSITE from the mapfile makes no difference.
>
> I strongly suspect I'm exporting the raster in the wrong type tiff,
> but it will take me a week to test them all and there could be an
> error elsewhere too.
>
> Can anyone see what I'm doing wrong?
>
> Z
>
>
> _______________________________________________
> grassuser mailing list
> grassuser at grass.itc.it
> http://grass.itc.it/mailman/listinfo/grassuser

I have had some similar problems wrt. exporting landsat data from GRASS -> 
mapsever. here is my approach:

#export tiled landsat-shade tiles:
for x in `seq 2 1 49`; do g.region rast=ca_landsat_shade_$x; r.out.tiff -t 
in=ca_landsat_shade_$x out=ca_landsat_shade_$x.tif; done

###BUG:
trying to add overviews to a tiff created with r.out.tiff:
gdaladdo ca_landsat_shade_20.tif.tif 2 4 8 16

0...10...20...30...40...50...60...70...80...90...100 - done.
ERROR 1: TIFFWriteEncodedTile:ca_landsat_shade_20.tif.tif: Tile 288 out of 
range, max 9
ERROR 1: TIFFWriteEncodedTile() failed.

# re-create geotiffs with gdal_translate (fix bug in tiffs created with 
r.out.tiff)
for x in `ls *.tif`; do gdal_translate -a_srs '+proj=aea +x_0=0.0 +y_0=0.0 
+lon_0=-96.0 +lat_0=40.0 +lat_1=20.0 +lat_2=60.0 +datum=NAD83' $x ../new/$x ; 
done 

#make overviews: (note that this will NOT work with images created from 
r.out.tiff)
for x in `ls *.tif`; do gdaladdo $x 2 4 8 16; done

terse i know, but hopefully helpful.

Dylan




-- 
Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
530.754.7341