Landsat image tiffs

[Zenon Panoussis]

Dylan Beaudette wrote:

> 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

You syntax confuses me and
http://mapserver.gis.umn.edu/docs/reference/mapfile/projection/
doesn't offer any help at all with it. What's x and y and why do you
have one lon and three lat? My map is 90N 90S 180W 180E wgs84.

> 2. valid projection string in your raster layer defs

As long as the tiff was UTM, putting the projection in the raster layer
didn't help. On the other hand, now that I warped the tiff to latlong,
its projection is the same as that of the rest of the map, so I shouldn't
need to declare it. When I try anyway with
    PROJECTION
      "proj=latlong"
      "ellps=WGS84"
      "datum=WGS84"
    END
in both the map and the layer, it makes no difference at all.

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

You mean that '[suffix]'? It's not there in the map, it was just a way
of referring in my posting to the three different tiffs I tried.

> 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

Hmm. grass says CELL for the raster, and I exporteded it with
'r.out.gdal type=UInt16'. This could theoretically be a source of errors,
but r.out.gdal can only export palettes in Byte and UInt16 mode (if I'm
not terribly mistaken) and r.out.tiff only exports 256 colours, so I should
use r.out.gdal if I want more than that. r.out.gdal converts CELL/Int32
to UInt16 and the output tiff says
  Band 1 Block=15706x1 Type=UInt16, ColorInterp=Palette
so I assume it is OK.

> 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

Nope, no gratification for me so far. The palette is there, has been there
all along, yet mapserver only displays a solid orange square:

Driver: GTiff/GeoTIFF
Size is 15706, 14422
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-77.361385,-10.633411)
Pixel Size = (0.00012981,-0.00012981)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( -77.3613848, -10.6334111) ( 77d21'40.99"W, 10d38'0.28"S)
Lower Left  ( -77.3613848, -12.5055129) ( 77d21'40.99"W, 12d30'19.85"S)
Upper Right ( -75.3226085, -10.6334111) ( 75d19'21.39"W, 10d38'0.28"S)
Lower Right ( -75.3226085, -12.5055129) ( 75d19'21.39"W, 12d30'19.85"S)
Center      ( -76.3419966, -11.5694620) ( 76d20'31.19"W, 11d34'10.06"S)
Band 1 Block=15706x1 Type=UInt16, ColorInterp=Palette
  Color Table (RGB with 65536 entries)
    0: 0,0,0,255
    1: 8,0,0,255
    2: 16,0,0,255
    3: 24,0,0,255
    4: 32,0,0,255
    5: 41,0,0,255
<snip 65530 lines>

I suspect that mapserver might be trying to interpret pixels instead of
just rendering them according to the image palette, but I don't know it
for sure, nor do I know how to change the behaviour.

BTW, what are those palette values? It looks like
<image pixel value>: R,G,B,<alpha?>
but it'd be nice to know what they stand for, not assume.

> 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

I'm not sure what this does. Did aea contain selected bands of one single
scene or did you wholesale import all the available bands?

> #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 

I didn't do that. I created the location by importing band 80, so as to
get a 14,25m base resolution, then imported bands 10, 20 and 30 in the
same location and run brovey on the lot. This should give me sharpened
RGB bands in 14,25 m resolution, and it seems it did:

r.info 007068.blue

 |   Type of Map:  raster               Number of Categories: 255
 |   Data Type:    DCELL
 |   Rows:         14461
 |   Columns:      15535
 |   Total Cells:  224651635
 |        Projection: UTM (zone 18)
 |            N: -1176408.75    S:   -1382478   Res: 14.25
 |            E:     464721    W:  243347.25   Res: 14.25
 |   Range of data:    min = 0.000000  max = 105.209790

You will notice that the data type is DCELL - wait a minute, I think
I found the pixel and skewing culprit. 14461x15535 it says here, one
too many, and the geometry is wrong. This gets fixed in the next step,
but this time I skipped i.landsat.rgb. last time I used it, so that's
pointing a finger in its direction (note: this refers to a related
discussion on the grass mailing list).

Anyway, I proceeded to merge $view.red $iew.green and $view.blue with
r.composite and got an output RGB raster like this:

 |   Data Type:    CELL
 |   Rows:         14460
 |   Columns:      15534
 |   Total Cells:  224621640
 |        Projection: UTM (zone 18)
 |            N: -1176415.875    S: -1382470.875   Res: 14.25
 |            E: 464713.875    W: 243354.375   Res: 14.25
 |   Range of data:    min = 0  max = 32767

with all the geometry back to what it's supposed to be and data type
CELL.

> 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 

Which i.landsat.rgb did you use (see
http://grass.itc.it/pipermail/grassuser/2006-August/035671.html and
the ensuing discussion)? The one in grass 6.1 doesn't take an 'out'
argument, it works on the input rasters (in my case the result of
brovey) and you have to run r.composite on them afterwards to get
your RGB final result.

> #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

I'm now working on one single scene, so this part is not yet relevant.
Later I plan landsat coverage over a bigger area, all of Peru, and I
fear that patched scenes will cause greater distorsions when reprojected
from UTM, than single scenes suffer individually. Anyway, that's a later
headache.

> ##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 

I gave up on r.out.tiff not only because of the bug, but also because
of its 8-bit limitation. r.out.gdal seems to do fine.

> #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'm skipping this until I've found a solution to the basic problem of
displaying the rasters in the first place.

> 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.

You'll have to hurry, or else I might beat you to it ;) Being new to GIS,
I find myself banging my head against the wall all the time, most of all
because of the tersity or absence of the documentation. That's particularly
true in grass: you get a maximally terse reference manual, not one single
example to go with it, you have to figure it all by trial and error and
endless googling. I am convinced that this oh so elusive error that has
taken me days and caused a long thread in two mailing lists is something
ridiculously simple, something silly that should be obvious, except it's
not obvious not to me. I'm sick of it, and the only thing I can do about
it is write some documentation myself as soon as I can do it decently.

Now the surprise: while I was typing that last paragraph, the problem got
solved. I have no idea why, but I do know how.

You will remember that I created tiffs both with r.out.tiff and r.out.gdal
and neither worked. I tried glal_translate to put back the geometry and
projection in the output of r.out.tiff, but it still didn't work. Mind
you, it was still in UTM.

The output of r.out.gdal had the geometry in it already. I tried warping
it to latlong and then it showed up on the map, but as a solid orange square
(note: square, its sides more or less straight west-east and north-south,
which is not the form of landsat scenes). The one thing I hadn't tried
was to correct the output of r.out.tiff with gdal_translate and then
gdalwarp the result of that. I realised this omission while answering your
posting and started working on it in the background. Just when I was ready
to hit send and go do something else to forget my frustrations, gdalwarp
finished and I tested the result. Bingo.

The correct procedure is the one I described in my opening posting, minus
i.landsat.rgb, up to the point of exporting. Then this:

r.out.tiff -t -l -v input=$view\_15me output=$view\_15meROT.tif

gdal_translate -a_srs '+proj=utm +zone=18 +datum=WGS84' \
    -a_ullr 243354.375 -1176415.875 464713.875 -1382470.875 \
    007068_15mROT.tif utm.tif

The values to -a_ullr are taken from the original NASA scene with

proj=utm
zone=`gdalinfo /vol/gis/sat/007068/p007r068_7p20000924_z18_nn80.tif \
   |grep PROJCS |awk '{print $6}' |sed 's/[^0-9]*//g'`
datum=WGS84
ul=`gdalinfo /vol/gis/sat/007068/p007r068_7p20000924_z18_nn80.tif \
   |grep 'Upper Left' |awk '{print $4}' |sed 's/,/ /' |sed 's/)//'`
lr=`gdalinfo /vol/gis/sat/007068/p007r068_7p20000924_z18_nn80.tif \
   |grep 'Lower Right' |awk '{print $4}' |sed 's/,/ /' |sed 's/)//'`
echo "gdal_translate -a_srs '+proj=$proj +zone=$zone +datum=$datum' \
   -a_ullr $ul$lr"

That's very ugly and could be much better, but it can serve as a basis
for later cosmetcs.

Having finished with gdal_translate, the last step is

gdalwarp -t_srs '+proj=latlong' utm.tif fixed.tif

The resulting tiff works with

  LAYER
    NAME "landsat"
    TYPE RASTER
    STATUS DEFAULT
    DATA "/vol/gis/lsat/fixed.tif"
#   MAXSCALE 100000
  END

in the map file and no PROJECTION statements anywhere.

It still remains to figure why the output of r.out.gdal doesn't work
the same way. Anyone with better understanding of tiff internals than
mine, please take a shot.

Z