0

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:

enter image description here

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.

6
  • what happen if you change dst_crs = rasterio.crs.CRS.from_dict(init='epsg:4326') to simply {'init': 'EPSG:4326'} ? Apr 1, 2021 at 8:19
  • @PierrickRambaud it still prints the destination array with 0 as shown in the post
    – ReutKeller
    Apr 1, 2021 at 8:24
  • When I plot it after the reprojection I get blank image
    – ReutKeller
    Apr 4, 2021 at 12:01
  • 1
    Why are you trying to reproject to the same CRS? Your question doesn't really make sense to me, what are you actually trying to do? Turn your numpy array into a rasterio dataset without writing to a file?
    – user2856
    Apr 5, 2021 at 1:49
  • @user2856 I have origianl raster, I read it as array, do some manipulations and then I want to clip it with shapefile I have. The problem is that the array loses the geographical data and then I can't clip it, unless I save it again as tiff and then open it again and clip it. I want to skip the step of save it as tiff and get the geographical data on my numpy array
    – ReutKeller
    Apr 5, 2021 at 9:14

2 Answers 2

3
+50

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'>
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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