Storing image with rasterio as gtiff modifies image histogram

Question

I'm trying to store and image as as gtiff, but the resulting image is changed, and I can't explain why.

I'm reading sentinel jp2 band images using rasterio, I rescale these images from [0-10.000] to [0-255]. Four bands are combined to create a [x,y,4] numpy array.

After that, the image is processed keeping the range [0-255] using numpy/opencv.

I want to store the resulting image in a gtiff file, but the image is modified during the process, and the original and resulting histograms are much different.

This is the profile from one of the original sentinel images:

{'driver': 'JP2OpenJPEG', 'dtype': 'uint16', 'nodata': None, 'width': 10980, 'height': 10980, 'count': 1, 'crs': CRS.from_dict(init='epsg:32630'), 'transform': Affine(10.0, 0.0, 300000.0, 0.0, -10.0, 4600020.0), 'blockxsize': 1024, 'blockysize': 1024, 'tiled': True}

To store the resulting image, I use this profile:

{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': 0, 'width': 10980, 'height': 10980, 'count': 4, 'crs': CRS.from_dict(init='epsg:32630'), 'transform': Affine(10.0, 0.0, 300000.0, 0.0, -10.0, 4600020.0), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'interleave': 'pixel', 'compress': 'jpeg'}

The code

# Here i'm reading one of the original sentinel images to get it's profile,  
# change it, and use it to geolocate the processed image
dataset = rasterio.open('one_of_the_original_sentinel_images')
profile = dataset.profile
dataset.close()

# change original profile to set geotiff driver
profile.update(
    **DefaultGTiffProfile(dtype=rasterio.uint8, count=img.shape[-1],
                          compress='jpeg', interleave='pixel',
                          # photometric='ycbcr', can't use it with 4-band image 
                          ))
# store the processed image (img) into a 4-band raster
with rasterio.open(img_filename, 'w', **profile) as dst:
    for ch in range(img.shape[-1]):
        # iterate over channels and write bands
        img_channel = img[:, :, ch]
        dst.write(img_channel, ch + 1)  # rasterio bands are 1-indexed

This is the image histogram before storing as gtiff with rasterio

Image histogram after storing it as gtif with rasterio (optained reading the gtif file directly with opencv.

UPDATE

Sorry, talking about the sentinel jp2, can be confusing mislead away from the main problem. I included it because I think the differences in profiles is the key of the problem.

My problem: I have an image stored in a numpy array, each channel in the range [0-255], and I want to geolocate it and store it in a gtiff file. To obtain he csr and affine, I use the original sentinel profile because I want the resulting image to have exactly the same geoposition.

The histograms I show above belong the image before and after storing it using rasterio. Besides, if I store the image directly using opencv, the image is right. So I think the pixel length or data range is not the cause of the problem.

So, which is the right way to get he width, height, affine and csr from a raster and use it to create a new raster with an image of the same size?

new_profile = DefaultGTiffProfile(dtype=rasterio.uint8, count=img.shape[-1], compress='jpeg',
                                      interleave='pixel')

# change the w/h, crs and affine from the sentinel original images
new_profile.update(crs=profile["crs"], transform=profile["transform"],
                   width=profile["width"],
                   height=profile["height"])

with rasterio.open(img_filename, 'w', **new_profile) as dst:
    for ch in range(img.shape[-1]):
        # iterate over channels and write bands
        img_channel = img[:, :, ch]
        dst.write(img_channel, ch + 1)  # rasterio bands are 1-indexed

Answer 1

You are changing many things (all of which will impact the histogram):

1. the pixel depth (16 bit to 8 bit),
2. Data range from 0-10 to 0-255,
3. file type from JP2 (potentially lossy wavelet compression) to Tiff (lossless since you are not using compression),
4. NoData from null to 0

You should, therefore, definitely expect the histogram to change because the output will certainly not be identical to the input. If you want the output to be identical - don't change it.

Answer 2

Using meta instead of profile solved the problem.

def geolocate_image(img, img_source, img_filename):
    """
    Creates a tiff raster to store the img param using geolocation
    information from another raster.

    :param img: Numpy array with shape [height, width, channels]
    :param img_source: raster to take the geolocation information from.
    :param img_filename: output raster filename
    :return:
    """

    with rasterio.Env():
        # read profile info from first file
        dataset = rasterio.open(img_source)
        meta = dataset.meta.copy()
        dataset.close()

        meta.update({"driver": "GTiff", "count": 4, 'dtype': 'uint8'})

        with rasterio.open(img_filename, 'w', **meta) as dst:
            for ch in range(img.shape[-1]):
                # iterate over channels and write bands
                img_channel = img[:, :, ch]
                dst.write(img_channel, ch + 1)  # rasterio bands are 1-indexed