1

I am trying to reproject a raster image from EPSG:4326 to EPSG:32737, and I believe I am having trouble because the original image is 16 bit- UINT16.

When I try and do the reprojection, I am ending up with an image that is all black, which suggests some possible truncation of the values at 255.

The code I am using is the same one from the Rasterio website. But I am not sure how to add the step of changing the datatype from Uint-16 to Uint-8 bit.

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'EPSG:32737'

with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    with rasterio.open('/tmp/RGB.byte.wgs84.tif', 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

I checked and the reproject function does not take the standard dtype argument. So I imagine this conversion has to happen in the src.meta metadata. But I was not sure if I can just change the setting, and that will take care of the conversion, or if I have to loop over the layers and change them, etc.

UPDATED:

So I wanted to show some of the output of this code. Mind that I tried this conversion with both a UINT8 image and a UINT16 image. The UINT8 image reprojection went just fine, but the UINT16 image had a problem. Here is an example:

This is a piece of the original raster drone image.

Original image before reprojection

And here is the image after using this same code:

Image after reprojection is all black

  • Are you (a) specifically trying to create a Byte output? (then rescale/stretch your values) or (b) unable to reproject a UInt16 and think that converting to Byte will help? (it won't) – user2856 Sep 26 at 5:19
  • I have these rasters that are about 38,000 x 38,000 pixels and I want to reproject to from 4326 to 32737. Now when I tried this the first time, using the code above, I got these really small files, like 21kb when the original file is 630mb. The code has worked before, and the only difference seems to be the UInt16 instead of UInt8. So I just need to get the reprojection to work. I suspect that UInt16 might mess things up further down the pipeline too, so was hoping to convert between data types. Do you have a better suggestion? – krishnab Sep 26 at 5:56
  • Integer rasters reproject with no issues, data type has nothing to do with it. I think you're conflating two issues/questions - "how do I reproject my raster?" and "how do I rescale my raster for final output". I'd focus on getting the reprojecting right (edit your question to remove Byte conversion stuff), then you can ask a follow up question about changing the datatype. – user2856 Sep 26 at 6:34
  • From your description I would say that you do not set the dst_crs right. You say that you want to warp to EPSG:32737 while you set your dst_crs to EPSG:4326 (see second to last line). – Johan Sep 26 at 6:52
  • @user2856 So I just tried this again. When I use the above code, I can reproject a raster that is UINT8 with no problem, but when I try it with a different raster that is UINT16, then I still end up with a black box. So barring any other reason, seems like data type is the cause. But I can perhaps update the question to show the resulting pictures of the black box versus original image. – krishnab Sep 26 at 7:01
1

I'm uncertain about rasterio but you can do this with gdal (which rasterio uses) in two steps.

1) use gdalwarp to warp from EPSG:4326 to EPSG:32737

2) use gdal_translate to scale the data from 16 bit to 8 bit.

With gdal, you can do the first step in memory (using vsimem)

In Python this would be:

import uuid
import gdal
import os

fn = "/path/to/some/file.tif"
outfn = "/path/to/some/outfile.tif"

xmin, ymin, xmax, ymax = # to be decided by the user
xres, yres = # to be decided by the user
extent = [xmin, ymin, xmax, ymax]
if not os.path.exists(outfn):

    # set the warp options
    warpoptions = dict(
        outputBounds=extent, # set the extent of the output
        xRes=xres, yRes=yres # set the x and y resolution of the output
        dstSRS="EPSG:32737", # destination SRS is EPSG:32737
        format="VRT" # warp to VRT
    )

    # set the translate options
    translateoptions = dict(
        scaleParams=[[0,2**16,0,2**8]], # rescale data from 16 to 8 bit
        outputType=gdal.GDT_Byte,
        creationOptions=["COMPRESS=LZW","TILED=YES"] # compress and tile the output
    )

    # step 1, warp input data to 32737 in memory
    desc = "/vsimem/" + uuid.uuid4().hex
    ds = gdal.Warp(desc, fn, **warpoptions)

    # step 2, translate the in memory dataset and apply the scaling
    ds2 = gdal.Translate(outfn, ds, **translateoptions)
    ds2 = None

    # close the in memory dataset
    ds = None
    gdal.Unlink(desc)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.