python - 在 geopandas 中移动阿拉斯加和夏威夷以获得 choropleths-6ren

python - 在 geopandas 中移动阿拉斯加和夏威夷以获得 choropleths

转载作者：行者123 更新时间：2023-12-05 02:00:00

在 R 中，我可以像这样移动阿拉斯加和夏威夷: https://www.storybench.org/how-to-shift-alaska-and-hawaii-below-the-lower-48-for-your-interactive-choropleth-map/

我可以像这样在 cartopy 中做到这一点: Showing Alaska and Hawaii in Cartopy map

我如何使用 geopandas 做类似的事情？

这是我所拥有的...

import requests
import geopandas
import numpy as np

np.random.seed(5)

URL = 'https://github.com/kjhealy/us-county/raw/master/data/geojson/gz_2010_us_040_00_500k.json'
us_json = requests.get(URL).json()

# save to disk as fname then load from file
usa = geopandas.read_file(fname)
usa["value"] = np.random.uniform(low=0,high=100, size=52)
usa.sort_values("NAME", inplace=True)
usa

现在绘制它:

usa.plot(column='value', cmap='Reds', scheme='quantiles', k=7, legend=True)

仅限大陆:

not_mainland = ['Alaska', 'Hawaii', "Puerto Rico"]
mainland_usa = usa.query('NAME not in @not_mainland')
mainland_usa.sort_values("NAME")
mainland_usa.plot(column='value', cmap='Reds', scheme='quantiles', k=7, legend=True)

我应该如何为阿拉斯加和夏威夷添加子图？

最佳答案

对于相邻的州，使用 geoDataFrame.plot() 照常进行。但阿拉斯加和夏威夷必须单独绘制为插图。在下面完整且可运行的代码中，有许多插入的注释解释了帮助读者的重要步骤。关于专题数据，我更喜欢使用更真实的人口数据(而不是问题中的妆容数据)来演示情节上的专题映射。包 mapclassify 用于对数据进行分类以进行适当的专题数据处理。

import matplotlib.pyplot as plt
import matplotlib
import cartopy
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import shapely.geometry as sgeom  #for box drawing
import geopandas as gpd
import numpy as np
import mapclassify as mc
import pandas as pd

import requests
import json

# population per sq-km
# include `Puerto Rico` but not used
popdensity = {
    'Alaska': 0.8,
    'District of Columbia': 4251,
    'Hawaii': 86,
    'Puerto Rico': 360,
    'New Jersey': 438.00,
    'Rhode Island': 387.35,
    'Massachusetts': 312.68,
    'Connecticut': 271.40,
    'Maryland': 209.23,
    'New York': 155.18,
    'Delaware': 154.87,
    'Florida': 114.43,
    'Ohio':  107.05,
    'Pennsylvania': 105.80,
    'Illinois': 86.27,
    'California': 83.85,
    'Virginia': 69.03,
    'Michigan': 67.55,
    'Indiana': 65.46,
    'North Carolina': 63.80,
    'Georgia': 54.59,
    'Tennessee': 53.29,
    'New Hampshire': 53.20,
    'South Carolina': 51.45,
    'Louisiana': 39.61,
    'Kentucky': 39.28,
    'Wisconsin': 38.13,
    'Washington': 34.20,
    'Alabama': 33.84,
    'Missouri': 31.36,
    'Texas': 30.75,
    'West Virginia': 29.00,
    'Vermont': 25.41,
    'Minnesota': 23.86,
    'Mississippi': 23.42,
    'Iowa': 20.22,
    'Arkansas': 19.82,
    'Oklahoma': 19.40,
    'Arizona': 17.43,
    'Colorado': 16.01,
    'Maine': 15.95,
    'Oregon': 13.76,
    'Kansas': 12.69,
    'Utah': 10.50,
    'Nebraska': 8.60,
    'Nevada': 7.03,
    'Idaho': 6.04,
    'New Mexico': 5.79,
    'South Dakota': 3.84,
    'North Dakota': 3.59,
    'Montana': 2.39,
    'Wyoming': 1.96}

# use this simple colormap
my_colormap = matplotlib.cm.Reds

# some settings
edgecolor = "gray"

# use this column for thematic mapping
theme_value = "pop_per_sqkm"

# A function that draws inset map
# ===============================
def add_insetmap(axes_extent, map_extent, state_name, facecolor, edgecolor, geometry):
    # create new axes, set its projection
    use_projection = ccrs.Mercator()      # preserves shape
    #use_projection = ccrs.PlateCarree()  # large distortion in E-W, bad for for Alaska
    geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
    sub_ax = plt.axes(axes_extent, projection=use_projection)  # normal units
    sub_ax.set_extent(map_extent, geodetic)  # map extents

    # option to add basic land, coastlines of the map
    # can comment out if you don't need them
    sub_ax.add_feature(cartopy.feature.LAND)
    sub_ax.coastlines()
    sub_ax.set_title(state_name)

    # add map `geometry`
    sub_ax.add_geometries([geometry], ccrs.PlateCarree(), \
                          facecolor=facecolor, edgecolor=edgecolor, lw=0.3)
    # +++ more features can be added here +++
    # plot box around the map
    extent_box = sgeom.box(map_extent[0], map_extent[2], map_extent[1], map_extent[3])
    sub_ax.add_geometries([extent_box], ccrs.PlateCarree(), color='none')

# access USA shapefile
# use data from internet
URL = 'https://github.com/kjhealy/us-county/raw/master/data/geojson/gz_2010_us_040_00_500k.json'

# more hard-code for settings
state_name = "NAME"  #defined column header
fname = "usa_json.json"

# Only on the first-run, set True in the next statement
if True:
    us_json = requests.get(URL).json()
    with open(fname, 'w') as file:
         file.write(json.dumps(us_json))

# saved it to disk as fname 
# then load from file
usa = gpd.read_file(fname)
usa.sort_values(state_name, inplace=True)

# make dataframe from population data `popdensity` dict object
usa_popden = pd.DataFrame.from_dict(popdensity, orient='index', 
                                  columns=['pop_per_sqkm'])
# merge `usa_popden` to main dataframe, `usa`
newusa = usa.merge(usa_popden, how='left', left_on='NAME', right_index=True)
# take only some columns in `newusa` for our operation
newusa = newusa[['NAME', 'pop_per_sqkm', 'geometry']]

# Data classification for thematic mapping
# choose number of classes of population density
# classes --> assigned colors in thematic mapping
num_classes = 7  #will get class-values: 0,1,2,3,4,5,6
sclass = mc.Quantiles(newusa["pop_per_sqkm"].values, k=num_classes)
#print(sclass)

# add new column, "sclass", for raw class values, sclass.yb
# its values will be used to assign color for polygon's facecolor
newusa["sclass"] = sclass.yb

# extract parts of the whole 'newusa' geodataframe for separate plotting/manipulation
# 'usa_main': excluding non-conterminous states
usa_main = newusa[~newusa[state_name].isin(["Alaska", "Hawaii", "Puerto Rico"])] # exclude these
#  re-project usa_main to equal-area conic projection "EPSG:2163"
usa_main.crs = {'init': 'epsg:4326'}
usa_main = usa_main.to_crs(epsg=2163)

# 'usa_more': non-conterminous states, namely, Alaska and Hawaii
usa_more = newusa[newusa[state_name].isin(["Alaska", "Hawaii"])]  # include these

# ------------ Plot --------------
# plot 1st part, using usa_main and grab its axis as 'ax2'
ax2 = usa_main.plot(column="sclass", legend=False, 
                    cmap=matplotlib.cm.Reds, ec=edgecolor, lw=0.4)

# manipulate colorbar/legend
fig = ax2.get_figure()
cax = fig.add_axes([0.9, .25, 0.02, 0.5])  #[left,bottom,width,height]
sm = plt.cm.ScalarMappable(cmap=my_colormap, 
        norm=plt.Normalize(vmin=min(newusa["sclass"]),vmax=max(newusa["sclass"])))

# clear the array of the scalar mappable
sm._A = []
cb = fig.colorbar(sm, cax=cax)
cb.set_label("Pop-density class")
# manipulate the axis seetings
ax2.set_frame_on(False)
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_xticklabels([])
ax2.set_yticklabels([])
ax2.set_title("Population Density Plot")

# add more features on ax2
# plot Alaska, Hawaii as inset maps
for index,state in usa_more.iterrows():

    if state[state_name] in ("Alaska", "Hawaii"):
        st_name = state[state_name]

        # set fill color, using normalized `sclass` on `my_colormap`
        facecolor = my_colormap( state["sclass"] / max(newusa["sclass"] ))

        if st_name == "Alaska":
            # (1) Alaska
            # Custom extent, relative size
            map_extent = (-178, -135, 46, 73)    # degrees: (lonmin,lonmax,latmin,latmax)
            axes_extent = (0.04, 0.06, 0.29, 0.275) # axes units: 0 to 1, (LLx,LLy,width,height)

        if st_name == "Hawaii":
            # (2) Hawaii
            # Custom extent, relative size
            map_extent = (-162, -152, 15, 25)
            axes_extent = (0.27, 0.06, 0.15, 0.15)

        # add inset maps
        add_insetmap(axes_extent, map_extent, st_name, \
                     facecolor, \
                     edgecolor, \
                     state["geometry"])

plt.show()

数据分类:

在上面的代码中，sclass 携带了人口密度的分类细节。要获取类列表，请运行以下代码:

low_val = 0
max_cval = len(sclass.bins)-1
print("Population density (per square_km)")
print("Class   Value_ranges")
for ix,val in enumerate(sclass.bins):
    #print(low_val, "< x <=", val)
    print("{ix:}    {low_val:.2f} < density <= {val:.2f}".format(ix=ix, low_val=low_val, val=val))
    #print(my_colormap(ix/max_cval))
    low_val = val

你应该得到下面的结果，上面的 map 应该伴随着它。

Population density (per square_km)
Class   Value_ranges
0    0.00 < density <= 7.48
1    7.48 < density <= 18.56
2    18.56 < density <= 30.50
3    30.50 < density <= 51.70
4    51.70 < density <= 75.38
5    75.38 < density <= 155.09
6    155.09 < density <= 4251.00

额外的颜色条操作

可以操纵颜色条将标签文本从简单的类编号 0,1,2,...6 更改为分类值范围。只需在上面代码的 plt.show() 行之前插入这些代码行，然后重新运行。

class_txts = []
low_val = 0
for ix,val in enumerate(sclass.bins):
    class_txts.append("{low_val:.1f}, {val:.1f}".format(ix=ix, low_val=low_val, val=val))
    low_val = val
cb.ax.set_yticklabels( class_txts )
cb.set_label("Pop-density range")

新情节将如下所示: