Landsat image tiffs

[Dylan Beaudette]

On 8/19/06, Zenon Panoussis <oracle at provocation.net> wrote:
> 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.
>

Well getting into GIS will take some reading. My projection happens to
be an Albers equal area projection, with 2 standard parallels. If
things do not make sense then you will need to get up to speed on the
conecpt of projection, datum, etc. This will have to be done outside
of the mapserver / GRASS documentation. If you are near a library, I
would recommend Principals of Cartography.


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

Here is what will make a difference:

1. your map file having a valid projection stanza
2. your raster layer having a valid projection stanza
3. you raster containing embedded projection data, or at least a world file

if any of these terms are new, please look them up.

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

ok.

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

I would ask on the GDAL mailing list about this. There has been a long
standing issue with the output from r.out.gdal.


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

I suspect an incorrect data type.

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

Looks like it.

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

this was just my import step - importing pre-warped images into GRASS
for color correction, mosaicing, and then re-tiling. I imported all
three bands for each scene in that step.

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

Ah. Well this is telling me that you arecreating DCELL maps- which is
not surprising coming from the fusion process. You will need to
quantize the images at some point back to CELL rasters.

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

I am not sure about your pixel shift-  was there was a mishap in the
interpretation of grid centered vs. pixel centered coordinates?

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

ok. there is the quantization that we need...

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

Ah. my notes are old. I was using the dev. version of that script.
glad to hear that the new one does not create new maps. please adjust
your commands as required.

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

Yes, you will have to deal with this at some point. do all of the
scenes come as UTM? where are you getting them from?


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

It is true that much of the documentation is terse. However, with a
solid understanding of the background material the documentation is
often sufficient. As with all FOSS projects, please feel free to
contribute to the documentation. I know that the GRASS team is always
looking for this type of help.

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

key concepts here, projection requires knowledge of sorce coordinate
system and destination. if either of these are not defined properly
you will not get the correct results.

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

Here is my interpretation. r.out.gdal is producing a geotiff encoded
in a way that mapserver cannot properly decode (orange rectangle). I
have found this to be the case for a while now. r.out.tiff will
produce acceptable output in 8-bit _or_ 24-bit color (note the -p
flag), however the file will be a TIFF and not a GeoTIFF.
gdal_translate can be used to add the spatial reference system, and
create a proper GeoTIFF form this output. This file can then be warped
and used in mapserver.

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

assign the output of each set of commands to a variable:

zone=`gdalinfo /vol/gis/sat/007068/p007r068_7p20000924_z18_nn80.tif \
 |grep PROJCS |awk '{print $6}' |sed 's/[^0-9]*//g'`

... +zone=$zone ...

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

Sure. the PROJECTION blocks are only needed if there is
1. a discrepancy in proejction between multiple layers
2. the projection cannot be ascertained from the input file of a layer

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

Dunno. I think that it is a bug with r.out.gdal. I would ask on the
GRASS list or GDAL or both.

Looking forward to seeing the working demo sometime!

Cheers,

dylan