1. Rasterio | Main
  2. Topics

Rasterio Python Library Merging R, G, B Bands to Create a GeoReferenced GeoTiff

nathan.raley@...
 

I am attempting to use rasterio to merge 3 bands from Sentinel-2 data, sourced as .JP2 files, into a georeferenced GeoTiff.  However, I am running into 2 problems doing this.

First of all, here is my code for the merge:
    #Setup the bands
    b4 = rio.open(path + Bands.B04.value[0])
    b3 = rio.open(path + Bands.B03.value[0])
    b2 = rio.open(path + Bands.B02.value[0])
   
    #Setup the path for the constructed RGB image
    tif = path + "RGB.tif"
   
    #Create the RGB image
    with rio.open(tif, 'w', driver='Gtiff', width=b4.width, height=b4.height, count=3, crs=b4.crs, tranform=b4.transform, dtype=b4.dtypes[0], photometric="RGB", tile=True, compress='lzw', bigtiff=True) as rgb:
        rgb.write(b2.read(1), 1)
        rgb.write(b3.read(1), 2)
        rgb.write(b4.read(1), 3)
        rgb.close()

Now, how do I copy the source projection information and bounds in the new GeoTiff I am creating?  Second, why is it showing up as all black when viewing the image in most common image browsers, but just fine in mapping software such as ArcMap or ArcGIS Pro?

Any ideas?
Sean Gillies
 

Hi,

On Wed, Mar 11, 2020 at 4:17 PM <nathan.raley@...> wrote:
I am attempting to use rasterio to merge 3 bands from Sentinel-2 data, sourced as .JP2 files, into a georeferenced GeoTiff.  However, I am running into 2 problems doing this.

First of all, here is my code for the merge:
    #Setup the bands
    b4 = rio.open(path + Bands.B04.value[0])
    b3 = rio.open(path + Bands.B03.value[0])
    b2 = rio.open(path + Bands.B02.value[0])
   
    #Setup the path for the constructed RGB image
    tif = path + "RGB.tif"
   
    #Create the RGB image
    with rio.open(tif, 'w', driver='Gtiff', width=b4.width, height=b4.height, count=3, crs=b4.crs, tranform=b4.transform, dtype=b4.dtypes[0], photometric="RGB", tile=True, compress='lzw', bigtiff=True) as rgb:
        rgb.write(b2.read(1), 1)
        rgb.write(b3.read(1), 2)
        rgb.write(b4.read(1), 3)
        rgb.close()
I have a couple comments on the code. You need to pass `tiled=True` to rasterio.open(), "tiled", and you should pass blockxsize and blockysize parameters as well. These must be multiples of 16. Also, you're passing "tranform" instead of "transform" and that will be ignored. This may explain why your output has no meaningful bounds.

I've found a pattern that makes dealing with a dataset's creation parameters a bit easier. In your case, it goes like this:

    rgb_profile = b4.profile
    rgb_profile.update(driver="GTiff", count=3, photometric="RGB", compress="LZW", bigtiff="IF_NEEDED", tiled=True, blockxsize=256, blockysize=256)
    with rasterio.open(tif, "w", **rgb_profile) as rgb:
        rgb.write(b2.read(1), 1)
        rgb.write(b3.read(1), 2)
        rgb.write(b4.read(1), 3)
 


Now, how do I copy the source projection information and bounds in the new GeoTiff I am creating?  Second, why is it showing up as all black when viewing the image in most common image browsers, but just fine in mapping software such as ArcMap or ArcGIS Pro?

Any ideas ?

I think it's likely that the common image browsers don't support some of the features of TIFF that GDAL uses. BIGTIFF is one of these.

--
Sean Gillies

1 - 2 of 2