在R中按阈值绘制人口栅格数据的解决方案咨询

Got it, I’ve run into this exact issue before—default raster plotting functions don’t let you easily split colors based on a hard threshold like population <1 vs >1. The fix involves creating a custom color mapping or layering your raster data to apply different colors to each subset. Below are two reliable approaches using Python tools that work seamlessly with GPW v4 data:

Approach 1: Custom Colormap with Boundary Normalization #

This method creates a single colormap that explicitly maps values ≤1 to red, and everything above to a gradient (we’ll use viridis here, but you can swap in any colormap you like).

Step-by-Step Code #

import rasterio
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
import numpy as np

# 1. Load your GPW raster file
# Replace the filename with your actual GPW .tif path
with rasterio.open('gpw_v4_population_count_rev10_2020_30_sec.tif') as src:
    pop_array = src.read(1)  # Read the first band (population count)
    raster_transform = src.transform  # Get spatial transform for plotting

# 2. Define your threshold and colors
population_threshold = 1
# Set red for low-population areas
red_color = "#ff0000"
# Grab a gradient colormap for high-population areas
high_pop_colors = plt.cm.viridis(np.linspace(0, 1, 256))
# Combine red with the gradient into a single custom colormap
custom_colors = np.vstack([np.array(plt.colors.to_rgb(red_color)), high_pop_colors])
custom_cmap = ListedColormap(custom_colors)

# 3. Create boundary normalization to split values at your threshold
bounds = [0, population_threshold, pop_array.max()]
norm = BoundaryNorm(bounds, custom_cmap.N)

# 4. Plot the raster
plt.figure(figsize=(12, 8))
im = plt.imshow(pop_array, cmap=custom_cmap, norm=norm, transform=raster_transform)

# Add a labeled colorbar
cbar = plt.colorbar(im, ticks=[0.5, (population_threshold + pop_array.max())/2])
cbar.ax.set_yticklabels(["Population ≤ 1", "Population > 1"])

plt.title("GPW Population Count: Red for Low-Density Areas")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

Approach 2: Layered Masked Plotting #

If you prefer more control (like using different color scales for each group), you can plot two masked versions of the raster: one for values ≤1 (red) and one for values >1 (your chosen gradient).

Step-by-Step Code #

import rasterio
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import numpy as np

# Load the raster data
with rasterio.open('gpw_v4_population_count_rev10_2020_30_sec.tif') as src:
    pop_array = src.read(1)
    raster_transform = src.transform
    # Handle NoData values (set to NaN so they don't plot)
    pop_array = np.where(pop_array == src.nodata, np.nan, pop_array)

# Create masks for each population group
low_pop_mask = pop_array <= 1
high_pop_mask = pop_array > 1

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

# First plot high-population areas with a gradient colormap
plt.imshow(
    np.where(high_pop_mask, pop_array, np.nan),
    cmap="viridis",
    transform=raster_transform
)

# Then overlay low-population areas in red
plt.imshow(
    np.where(low_pop_mask, 1, np.nan),
    cmap=ListedColormap(["#ff0000"]),
    transform=raster_transform
)

# Add colorbar for high-population values
cbar = plt.colorbar(plt.cm.ScalarMappable(cmap="viridis"), ax=plt.gca())
cbar.set_label("Population > 1")

# Add a legend for the red low-population areas
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor="#ff0000", label="Population ≤ 1")]
plt.legend(handles=legend_elements, loc="upper right")

plt.title("GPW Population Count: Layered Threshold Coloring")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

Key Notes #

• NoData Handling: Both examples account for NoData values (common in GPW data) by converting them to np.nan so they don’t interfere with your color mapping.
• Customization: Swap out "viridis" for any other Matplotlib colormap (like "plasma" or "magma") to change the high-population color scheme. You can also adjust the red hex code to any shade you prefer.
• Why Default plot() Fails: The built-in plot() functions for rasters use continuous color scales that map every value to a color in sequence. They don’t support hard threshold splits out of the box, which is why we need to build custom mappings or layered plots.

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