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
