On 8/19/06, Zenon Panoussis <[EMAIL PROTECTED]> 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
