Landsat image tiffs

[Dylan Beaudette]

On Friday 18 August 2006 23:43, Zenon Panoussis wrote:
> Zenon Panoussis wrote:
> > In a map that displays etopo and srtm tiffs without problems,
> > I have this layer:
> >
> >  LAYER
> >    NAME "landsat"
> >    TYPE RASTER
> >    STATUS DEFAULT
> >    DATA "007068_15me[suffix].tif"
> >  END
> >
> > mapserver refuses to display this particular layer no matter what I do.
>
> The base map uses UNITS DD. Experimenting on, I added
>
> PROJECTION
>   "proj=utm"
>   "ellps=WGS84"
>   "datum=WGS84"
>   "zone=18"
>   "units=m"
>   "south"
>   "no_defs"
> END
>
> to the layer. The layer was still not displayed. I then went in the
> opposite direction, removed the projection definition from the layer
> and run
>
>   gdalwarp -t_srs '+proj=latlong' old.tif new.tif
>
> Now mapserver actually displays the raster, but only as a solid-colour
> orange square.
>
> I've already been barking up so many trees trying to get this to work,
> and there's still a whole forest left. Does anyone have a sample layer
> and creation command for landsat images, and would be willing to post
> it here?
>
> Z

Zenon,

I would double check a couple of things. 

1. valid projection settings for your MAP file:

#my sample:
PROJECTION
"+proj=aea +x_0=0.0 +y_0=0 +lon_0=-96 +lat_0=40.0 +lat_1=20 +lat_2=60.0 
+datum=NAD83"
END

2. valid projection string in your raster layer defs

3. remove mapserver variable  from data source name:
 DATA "007068_15me[suffix].tif"
                                   ^^^^^^^

4. verify your datatypes in GRASS when creating / compositing mulitple rasters 
i.e. use r.info to check CELL, FCELL, DCELL
note that these will be equiv. to Int32, Float32, and Float64 in GDAL land

5. try and output an 'image' with palette, rather than the raw data with 
r.out.tiff - by having the color table built-in to the tiff, you will get 
instant gratification in mapserver

here is a list of how i previously processed landsat data: note that the 
scenes were imported as separate bands, resulting in 3 CELL type raster files 
for each scene.

#import landsat scenes
cd aea
for x in `ls *.tif`; do r.in.gdal in=$x out=`basename $x .tif`; done

#make composite landsat scenes from 3 bands:
#recall that this will be done in the native resolution of each set of images
#note special syntax to remove the tailing character from the 
for x in `g.mlist rast pattern="*_nn1"`
do 
g.region rast=$x 
echo "working on image: ${x%?}"
#r.composite r=${x%?}3 g=${x%?}2 b=${x%?}1 out=${x%?}.rgb 
#new better method
i.landsat.rgb red=${x%?}3 green=${x%?}2 blue=${x%?}1 out=${x%?}.rgb 
done 

#set 0 = NULL in the composite landsat scenes
for x in `g.mlist rast pattern="*.rgb"`; do g.region rast=$x ; r.null map=$x 
setnull=0; done

#zoom to the full extent
g.region -a rast=az_30m_dem res=30

#patch together the landsat scenes
r.patch in=`g.mlist type=rast pattern=*.rgb sep=","` out=az_landsat

...

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

# re-create geotiffs with gdal_translate (fix bug in tiffs created with 
r.out.tiff)
mkdir new

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)
cd new
for x in `ls *.tif`; do gdaladdo $x 2 4 8 16; done


I will put a more detailed version of these general steps online for others 
who might be interested in mass blending / tiling / export to mapserver.

Cheers,

Dylan



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