Combining Landsat Bands

[AtmoPierce]

Hello,

This may be a general GDAL question, but I'll give it a shot here.

My initial goal is to be able to combine RGB bands from a Landsat 9 image. Here's is my generic test code. Having general issues with the striping, and color.

`
let red_dataset = Dataset::open(red_path).unwrap();
let green_dataset = Dataset::open(green_path).unwrap();
let blue_dataset = Dataset::open(blue_path).unwrap();

let red_rasterband: RasterBand = red_dataset.rasterband(1).unwrap();    
let green_rasterband: RasterBand = green_dataset.rasterband(1).unwrap();
let blue_rasterband: RasterBand = blue_dataset.rasterband(1).unwrap();

let red_block_size = red_rasterband.block_size();
let green_block_size = green_rasterband.block_size();
let blue_block_size = blue_rasterband.block_size();

println!("Block Size: {:?}", red_block_size);
println!("Band Interpretation: {:?}", red_rasterband.color_interpretation());
println!("Red Max/Min {:?}", red_rasterband.compute_raster_min_max(true).unwrap());

let image_size = red_dataset.raster_size();

let image_size_u32: (u32, u32) = (image_size.0.try_into().unwrap(), image_size.1.try_into().unwrap());
let mut img = RgbImage::new(image_size_u32.0, image_size_u32.1); 

if let Ok(red_rv) = red_rasterband.read_as::<u8>((0, 0),            image_size, red_block_size,  Some(ResampleAlg::CubicSpline)) {
    if let Ok(green_rv) = green_rasterband.read_as::<u8>((0, 0),    image_size, green_block_size,Some(ResampleAlg::CubicSpline)) {
        if let Ok(blue_rv) = blue_rasterband.read_as::<u8>((0, 0),  image_size, blue_block_size, Some(ResampleAlg::CubicSpline)) {
            for pixel_place_x in 0..image_size_u32.0{
                for pixel_place_y in 0..image_size_u32.1{
                    let access_point: usize = (pixel_place_x + pixel_place_y).try_into().unwrap();
                    img.put_pixel(pixel_place_x, pixel_place_y, Rgb([red_rv.data[access_point], green_rv.data[access_point], blue_rv.data[access_point]]));
                }
            }
            //println!("{:?}", rv.data);
            img.save("results/test.png").unwrap();
        }
    }
}`

The output file:

My thoughts, something going on with the block size, something going on with my access point to the strided vector, or some kind of scaling in the TIFF. Does anybody have any suggestions?

[lnicola]

It's because of your index computation. The way you're doing it now, pixels (0, 10) and (10, 0) will map to the same access_point, 10. This causes the diagonal lines -- the diagonals have the same values because they're actually filled with the same pixel.

So here you need to think about the pixel order. Since they're laid out in rows (first row, then the second row etc.), pixel (x, y) will be at y * width + x. Note, though, that some language and libraries will store 2D arrays by columns instead of rows.

Finally, it won't affect the result but, unless the compiler is being very smart, it will be faster if you access the pixels in a sequential pattern. You're looping by x, then by y, which means you're iterating odero columns. You should switch the two for loops, so that the x varies faster.

[AtmoPierce]

That worked great. Shall be working on colors next.

Results:

[lnicola]

That's because the values in that L9 raster are larger than 255 (the u8 limit). You can download e.g. QGIS and use the Identify tool to check the pixel values. You'll want at least to choose an input range (e.g. 0..10000), then map it down to 0..255. It's somewhat like Levels in any photo editor (except there's no midpoint slider). Once you get that working, you can try applying gamma correction, since the L9 values are linear and the images will appear too dark otherwise.

[AtmoPierce]

Mean for this image is at like 13000, will definitely be doing some mapping to fix contrast/brightness. I was able to implement something similar in python before, but it definitely used some high-level packages. Will be looking forward to implementing this! Since I have access to the raw TIFF values should be able to directly manipulate them. Here's a quick and dirty spline (wouldn't call this real data):