Why the extents of GeoDataFrame and masked file are different?

Question

vector = some shape file, MultiPolygon in GeoDataFrame. File crs and vector crs are the same.

vector.geometry

0 MULTIPOLYGON (((47.59780 -19.38530, 47.59330 -... Name: geometry, dtype: geometry

extent = [vector.total_bounds[0], vector.total_bounds[2], vector.total_bounds[1], vector.total_bounds[3]]

with rio.open(URL) as file:
    crs = vector.crs
    arr, out_transform = mask(file, vector.geometry, crop=True)

shapes = rio.features.shapes(arr, transform=out_transform)
geometry = []
vals = []
for polygon, value in shapes:
    if value > 0:
        vals.append(value)
        geometry.append(shape(polygon))

gdf2 = gpd.GeoDataFrame({'vals': vals, 'geometry': geometry}, crs=crs)

Now, when I plot the masked raster and resulting GeoDataFrame, I expect grids to match each other, but the reality is somewhat different...

fig, ax = plt.subplots(figsize=(8,8))
arr[arr <= 0 ] = 'nan'
base = show(arr[0], extent=extent, ax=ax)
colorbar = base.get_images()[0]
fig.colorbar(colorbar, ax=ax)
vector.plot(facecolor='none', edgecolor='grey', ax=ax)

gdf2.iloc[:].plot(facecolor='none', edgecolor='red', ax=ax);

I find, that the extents of GeoDataFrames are different for some reasons.

print(vector.total_bounds)
print(gdf2.total_bounds)

[ 46.495 -19.5584 48.037 -17.707 ] [ 46.744999 -19.540326 48.078327 -17.873666]

And this is despite the specified transform parameter: out_transform

Affine(0.083333, 0.0, 46.415761, 0.0, -0.083333, -17.666236000000012)

What do I do wrong?

Libraries imported:

import numpy as np
import rasterio as rio
from rasterio.mask import mask
from rasterio.plot import show
from shapely.geometry import shape
import matplotlib.pyplot as plt
import geopandas as gpd

Answer 1

You need to pass a transform to rasterio.plot.show:

base = show(arr[0], extent=extent, ax=ax, transform=out_transform)

Reproducible example:

import numpy as np
import rasterio as rio
from rasterio.mask import mask
from rasterio.plot import show
from shapely.geometry import shape
import matplotlib.pyplot as plt
import geopandas as gpd


# Some dummy data
url = "https://github.com/rasterio/rasterio/raw/main/tests/data/RGB2.byte.tif"
features = {
    "type": "FeatureCollection",
    "features": [
        {
            "id": "0",
            "type": "Feature",
            "properties": {"id": 1},
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [[238793, 2653121], [243276, 2656230], [249353, 2651980], [242146, 2649325], [238793, 2653121]]]
            }
        }
    ]
}

with rio.open(url) as raster:
    crs = raster.crs

    vector = gpd.GeoDataFrame.from_features(features, crs)
    extent = [vector.total_bounds[0], vector.total_bounds[2], vector.total_bounds[1], vector.total_bounds[3]]

    arr, out_transform = mask(raster, vector.geometry, crop=True)

shapes = rio.features.shapes(arr, transform=out_transform)
geometry = []
vals = []
for polygon, value in shapes:
    if value > 0:
        vals.append(value)
        geometry.append(shape(polygon))

gdf2 = gpd.GeoDataFrame({'vals': vals, 'geometry': geometry}, crs=crs)

fig, ax = plt.subplots(figsize=(8,8))
arr[arr <= 0 ] = 0  # replaced nan with 0 just to get the code
                    # to work with above byte data
base = show(arr[0], extent=extent, ax=ax, transform=out_transform)
colorbar = base.get_images()[0]
fig.colorbar(colorbar, ax=ax)
vector.plot(facecolor='none', edgecolor='grey', ax=ax)

gdf2.iloc[:].plot(facecolor='none', edgecolor='red', ax=ax);

plt.show()