1

I have a shapefile that I use to crop GeoTIFFs. However everytime I use my code to do that, it slightly shifts my entire cropped TIFF to the East and I have no idea why. The pictures below show how well the original TIFF (uncropped) fits, and the second image shows how shifted the cropped TIFF is (I plotted only the boundaries of my shapefile, so it is more readable). TIFFs and shapefile are all in EPSG:32607 (UTM 7N).

Link to a Google Drive containing: .shp file, .shx file, and the TIFF.

Google_drive_link

Do you know why that is happening ?

Original uncropped tif Cropped/shifted tif

The code that I use is:

def Clipper(raster, option):

  path = './glacier.shp'
  
  # Read Shapefile
  with fiona.open(path, "r") as shapefile:
      shapes = [feature["geometry"] for feature in shapefile]

  # read imagery file
  out_image, out_transform = rio.mask.mask(raster, shapes, crop=True, nodata = np.nan)

  # Check that after the clip, the image is not empty
  test = out_image[~np.isnan(out_image)]

  if test[test>0].shape[0] == 0:
    passaran = 0
  else:
    passaran = 1

  # convert the zeros into NaNs
  
  return(out_image, passaran)

# Read a dummy src
raster = rio.open('./S2_2018_09_08.tif')
out_meta = raster.meta

# Clip the raster
array_out, passaran = Clipper(raster)

# Save the clip as a tif
with rio.open('./clipped.tif', "w", **out_meta) as dest:
              dest.write(array_out)

Metadata raster: enter image description here Metadata shapefile: enter image description here

4
  • Try using the array shape and out_transform as per the docs example instead of raster = rio.open('/content/dummy.tif')
    – user2856
    Apr 23, 2022 at 0:15
  • It doesn't work, it actually shifts it even more. The only option that shifts it the least is what I do, without updating the metadata
    – Nihilum
    Apr 23, 2022 at 0:33
  • I'm so sorry, I added a link to a Google Drive containing the files. Thank you for your time !!
    – Nihilum
    Apr 23, 2022 at 3:39
  • Thank you for trying ! I wonder if it comes from my installation. I work with Google Colab so I have no idea if it somehow impacts things...
    – Nihilum
    Apr 23, 2022 at 4:34

1 Answer 1

3

Your code worked for me with minor modifications:

import rasterio as rio, fiona
import rasterio.mask as mask
import numpy as np

raster = 'S2_2018_09_08.tif'
clipped = 'clipped.tif'
path = 'glacier.shp'


def Clipper(raster, vector):


    # Read Shapefile
    with fiona.open(vector, "r") as shapefile:
        shapes = [feature["geometry"] for feature in shapefile]

    with rio.open(raster) as src:
        # read imagery file
        out_image, out_transform = mask.mask(src, shapes, crop=True, nodata=np.nan)

        # Check that after the clip, the image is not empty
        test = out_image[~np.isnan(out_image)]

        if test[test > 0].shape[0] == 0:
            raise RuntimeError("Empty output")

        out_meta = src.profile
        out_meta.update({"height": out_image.shape[1],
                         "width": out_image.shape[2],
                         "transform": out_transform})

    return (out_image, out_meta)


# Clip the raster
array_out, out_meta = Clipper(raster, path)

# Save the clip as a tif
with rio.open(clipped, "w", **out_meta) as dest:
    dest.write(array_out)

enter image description here

enter image description here

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.