如何在Python中将自定义TIFF图像作为月球南极网格的底图并叠加多幅TIFF图像

First off, let’s fix a quick issue in your original code—you’re using npstere (north polar stereographic) but need spstere for the south pole. Then we’ll walk through adding support for your custom TIFF base map and overlays.

Step 1: Install Required Libraries #

You’ll need rasterio to handle geospatial TIFF files properly. Install it alongside your existing dependencies with:

pip install rasterio matplotlib basemap numpy

Step 2: Full Working Code #

Here’s the corrected, expanded code that loads your base TIFF and lets you overlay additional TIFFs with transparency:

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.warp import reproject, Resampling
from rasterio.transform import from_bounds

# Setup South Polar Stereographic Basemap
# Use 'spstere' for south pole; boundinglat sets the visible edge (e.g., -60°S for wider view)
m = Basemap(projection='spstere', boundinglat=-60, lon_0=0, resolution='l')

def plot_tiff_on_basemap(m, tiff_path, alpha=1.0):
    """Helper to load and plot TIFFs on the Basemap, handling coordinate conversion/reprojection."""
    with rasterio.open(tiff_path) as src:
        # Read first band (adjust to [1,2,3] if using RGB TIFFs)
        tiff_data = src.read(1)
        src_crs = src.crs
        src_bounds = src.bounds

        # If TIFF uses lon/lat coordinates, convert to Basemap's projected system
        if src_crs.is_geographic:
            x_left, y_bottom = m(src_bounds.left, src_bounds.bottom)
            x_right, y_top = m(src_bounds.right, src_bounds.top)
            m.imshow(tiff_data, extent=[x_left, x_right, y_bottom, y_top], origin='upper', alpha=alpha)
        else:
            # Reproject TIFF to match Basemap's South Polar Stereographic projection (EPSG:3031)
            target_crs = 'EPSG:3031'
            if src_crs != target_crs:
                tiff_data, _ = reproject(
                    tiff_data,
                    src_transform=src.transform,
                    src_crs=src_crs,
                    dst_transform=from_bounds(*m.bounds, width=tiff_data.shape[1], height=tiff_data.shape[0]),
                    dst_crs=target_crs,
                    resampling=Resampling.bilinear
                )
            m.imshow(tiff_data, extent=m.bounds, origin='upper', alpha=alpha)

# Plot your base map TIFF
plot_tiff_on_basemap(m, 'your_base_map.tif')

# Overlay additional TIFFs (adjust alpha for transparency)
plot_tiff_on_basemap(m, 'overlay1.tif', alpha=0.5)
# Add more overlays as needed:
# plot_tiff_on_basemap(m, 'overlay2.tif', alpha=0.3)

# Draw grid lines tailored to the south pole
m.drawparallels(np.arange(-90., 0., 20.), labels=[1,0,0,0])
m.drawmeridians(np.arange(-180., 180., 20.), labels=[0,0,0,1])

plt.title("Lunar South Polar Grid with Custom Base Map & Overlays")
plt.show()

Key Details to Note: #

• Projection Fix: Switched npstere to spstere to correctly focus on the south pole. Tweak boundinglat (e.g., to -80) for a tighter view of the lunar pole.
• TIFF Flexibility: The helper function works with both lon/lat and projected TIFFs, automatically converting coordinates or reprojecting to match the Basemap.
• Transparency: Use the alpha parameter when calling plot_tiff_on_basemap to make overlays semi-transparent, so your base map stays visible underneath.
• RGB TIFFs: If your base map is RGB, modify src.read(1) to src.read([1,2,3]) and the imshow call will handle the color data.

Alternative: Cartopy (Recommended for Long-Term Use) #

Since Basemap is deprecated, consider switching to Cartopy for better ongoing support. Here’s a quick snippet to get you started:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import rasterio

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.SouthPolarStereo()})
ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())

# Load and plot base map
with rasterio.open('your_base_map.tif') as src:
    data = src.read()
    ax.imshow(data[0], extent=src.bounds, transform=ccrs.PlateCarree(), origin='upper')

# Add overlays using the same pattern
plt.show()

内容的提问来源于stack exchange，提问作者Nidhi Verma