Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Geotiff that is located out of bounds of the specified coordinate system is incorrectly resized #2320

Open
ilopata1 opened this issue Jan 29, 2024 · 0 comments

Comments

@ilopata1
Copy link

ilopata1 commented Jan 29, 2024

Description

I have a Geotiff with a CRS of NAD83 / Wisconsin South (ftUS). The file can be downloaded from
this link.

However, the CRS was incorrectly assigned to this image. The actual location was in NAD83 / Illinois West (ftUS).

When I attempt to read and display this image using the code below, the image is incorrectly resized and, as a consequence, the image is not displayed (view limits and data limits are consequentially miscalculated also)

I understand that the correct CRS should have been assigned. If it had been, the image would, of course, display correctly. However, I would expect a warning of some kind instead of a failure with this image.

Code to reproduce

import os
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import rasterio
import numpy as np

fname = r'image.tif'

fig = plt.figure(figsize=(8, 12))

with rasterio.Env():
    with rasterio.open(fname, 'r') as src:

        #Get the epsg code from the source file
        code = src.crs.to_epsg()
       
        try:
            projection = ccrs.epsg(code)
        except:
            raise ValueError("Error converting EPSG code")

        img = src.read()
        img = np.transpose(img, [1,2,0])

        ax = plt.axes(projection=projection)

        # Fetch extent of raster in projection units
        xmin, ymin, xmax, ymax = src.bounds
        print("Image size in data coords before imshow",xmax - xmin,"x",ymax - ymin)
        img_extent = [xmin,xmax,ymin,ymax]
        
        # Add the image.     
        image = ax.imshow(img, origin='upper', extent=img_extent, transform=projection)

        xmin, xmax, ymin, ymax = image.get_extent()
        print("Image size in data coords after imshow", xmax - xmin,"x",ymax - ymin)
        
        plt.show()

I believe this may be a bug. If I apply the following change to geoaxes.py the issue is resolved, but I have no idea what other problems such a change may cause.

index 18e38a4..13643d8 100644
--- "a/geoaxes.py"
+++ "b/geoaxes.py"
@@ -1310,7 +1310,8 @@ class GeoAxes(matplotlib.axes.Axes):
                 raise ValueError('Expected a projection subclass. Cannot '
                                  'handle a %s in imshow.' % type(transform))

-            target_extent = self.get_extent(self.projection)
+            target_extent = img.get_extent(self.projection)
             regrid_shape = kwargs.pop('regrid_shape', 750)
             regrid_shape = self._regrid_shape_aspect(regrid_shape,
                                                      target_extent)
Full environment definition

Operating system

Windows 11

Cartopy version

Python=3.9
Cartopy=0.22
Matplotlib=3.8.2
PyProj=3.6.1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant