Set CRS to numpy array with no geographical data using rasterio

Question

I have 1 dimension numpy array that has no geographical information. The numpy array was originally satellite image that I have manipulated in order to get this 1 dimension array.
I have the crs and transform of this array.
However, until now, in order to have geographic data on this array I had to write it as tiff and only then I could get back the geographical data and do things such as crop the raster and more.

I know there is option to reproject an array with no geographical data, as mentioned here . I have tried to do this with my data as following, based on the doc:

src_shape = (886, 1138)
rows, cols = (886, 1138)
d = 1.0/240 # decimal degrees per pixel
# The following is equivalent to
# A(d, 0, -cols*d/2, 0, -d, rows*d/2).
src_transform = A.translation(-cols*d/2, rows*d/2) * A.scale(d, -d)
src_crs = {'init': 'EPSG:4326'}
source = img_pred
dst_shape = (886, 1138)
dst_transform = rasterio.transform.from_bounds(*bbox, width= bbox_size[1], height= bbox_size[0])
this gets the transform which is in this case : Affine(0.00011618510158013974, 0.0, -44.75374, 0.0, -7.050087873462236e-05, -6.80845)
dst_crs = rasterio.crs.CRS.from_dict(init='epsg:4326')
destination = np.zeros(dst_shape, np.uint8)
reproject(
    img_pred,
    destination,
    src_transform=src_transform,
    src_crs=src_crs,
    dst_transform=dst_transform,
    dst_crs=dst_crs,
    resampling=Resampling.nearest)

Then I get new array with values 0:

So my question is : how do I get my array with the values reprojected correctly?

I assume that the problem might be in the definition of the original CRS, (src_crs), as my numpy array has no crs at all, but not sure if this is the reason.

Answer 1

Reprojecting is not the way. You need to create a rasterio dataset. You can create an in-memory dataset using a rasterio.MemoryFile using the georeferencing from your original dataset.

Assuming you have not clipped or resampled your numpy array, here is a method to generate an in-memory dataset:

from contextlib import contextmanager
import rasterio
use context manager so DatasetReader and MemoryFile get cleaned up automatically
@contextmanager
def mem_raster(array, source_ds):
    profile = source_ds.profile
    profile.update(driver='GTiff', dtype=array.dtype)
with rasterio.MemoryFile() as memfile:
    with memfile.open(**profile) as dataset: # Open as DatasetWriter
        dataset.write(array)

    with memfile.open() as dataset:  # Reopen as DatasetReader
        yield dataset  # Note yield not return



with rasterio.open('path/to/raster') as src:
    array = src.read() / 2.0
    with mem_raster(array, src) as mem:
        print(mem.dtypes, mem.crs, mem.bounds)
        print(repr(mem))

Which for my data outputs:

('float32',) EPSG:4326 BoundingBox(left=147.2, bottom=-35.54, right=147.7, top=-34.54)
<open DatasetReader name='/vsimem/93fd15a3-727d-4186-aeed-553a1cda1899/93fd15a3-727d-4186-aeed-553a1cda1899.tif' mode='r'>

Answer 2

This Page seems to have the answer you are looking for. You need to use rasterio.warp.calculate_default_transform()

1. Create your raster in EPSG:4326 as you have already done

2. reproject according to the method in the link (just tested it - it works)

    from rasterio.warp import calculate_default_transform, reproject, Resampling
   filename = 'my_4326_file.tif'
    with rasterio.open(filename) 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
        })
        outfilename = '/Users/zoid/Desktop/my_32755_file.tif'
    with rasterio.open(outfilename, '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.lanczos)
   
   

This resulted in a transformed raster that matches the untransformed one.