0

我按照亚当·赛明顿（Adam Symington）的这个优秀指南，成功创建了以下沙巴（马来西亚的一个州，东南亚国家）的地形图。左上角尴尬的黑色斑点是我试图在地图上绘制某些坐标的尝试。

在此处输入图像描述

我想通过以下方式改进此图表：

编辑：我已经弄清楚了第（1）项并在下面发布了解决方案。(2) 和 (3) 待定。

  1. [已解决]数据框包含 该sch州所有学校的坐标。我想在地图上标出这些。我怀疑它目前正在变得不稳定，因为轴不是“地理轴”（意思是，不使用纬度/经度比例） - 您可以通过设置来确认这一点ax.axis('on')。我该如何解决这个问题？ [解决了]

  2. 我想将实际领土之外的部分设置为白色。通话ax.set_facecolor('white')不起作用。我知道将其设置为灰色的具体内容是ax.imshow(hillshade, cmap='Greys', alpha=0.3)线条（因为更改 cmap 会更改背景）；我只是不知道如何改变它，同时保持地图中的颜色为灰色。

  3. 如果可能的话，我希望地图的轮廓是黑色的，但这只是迂腐。

重现上图的所有代码如下。该downloadSrc函数在本地文件夹中获取并保存依赖项（包含地形数据的 5.7MB 二进制文件和包含要绘制的点坐标的 0.05MB csv）；你只需要运行一次。

import rasterio
from rasterio import mask as msk

import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.colors import ListedColormap

import numpy as np
import pandas as pd
import geopandas as gpd
import earthpy.spatial as es
from shapely.geometry import Point

def downloadSrc(dl=1):
    if dl == 1:
        import os
        os.mkdir('sabah')

        import requests
        r = requests.get('https://raw.githubusercontent.com/Thevesh/Display/master/sabah_tiff.npy')
        with open('sabah/sabah_topog.npy', 'wb') as f: f.write(r.content)

        df = pd.read_csv('https://raw.githubusercontent.com/Thevesh/Display/master/schools.csv')
        df.to_csv('sabah/sabah_schools.csv')

# Set dl = 0 after first run; the files will be in your current working directory + /sabah
downloadSrc(dl=1)

# Load topography of Sabah, pre-saved from clipped tiff file (as per Adam Symington guide)
value_range = 4049
sabah_topography = np.load('sabah/sabah_topog.npy')

# Load coordinates of schools in Sabah
crs={'init':'epsg:4326'}
sch = pd.read_csv('sabah/sabah_schools.csv',usecols=['lat','lon'])
geometry = [Point(xy) for xy in zip(sch.lon, sch.lat)]
schools = gpd.GeoDataFrame(sch, crs=crs, geometry=geometry)

# Replicated directly from guide, with own modifications only to colours
sabah_colormap = LinearSegmentedColormap.from_list('sabah', ['lightgray', '#e6757b', '#CD212A', '#CD212A'], N=value_range)
background_color = np.array([1,1,1,1])
newcolors = sabah_colormap(np.linspace(0, 1, value_range))
newcolors = np.vstack((newcolors, background_color))
sabah_colormap = ListedColormap(newcolors)
hillshade = es.hillshade(sabah_topography[0], azimuth=180, altitude=1)

# Plot
plt.rcParams["figure.figsize"] = [5,5]
plt.rcParams["figure.autolayout"] = True
fig, ax = plt.subplots()
ax.imshow(sabah_topography[0], cmap=sabah_colormap)
ax.imshow(hillshade, cmap='Greys', alpha=0.3)
schools.plot(color='black', marker='x', markersize=10,ax=ax)

ax.axis('off')
plt.show()
4

1 回答 1

0

事实证明，我已经给了自己回答点（1）的提示，并且还设法解决了（2）。

对于 (1)，点只需要重新缩放，我们得到：

在此处输入图像描述

我通过从底层 shapefile 获取地图的最大/最小点，然后根据轴的最大/最小点对其进行缩放，如下所示：

# Get limit points
l = gpd.read_file('param_geo/sabah.shp')['geometry'].bounds
lat_min,lat_max,lon_min,lon_max  = l['miny'].iloc[0], l['maxy'].iloc[0], l['minx'].iloc[0], l['maxx'].iloc[0]
xmin,xmax = ax.get_xlim()
ymin,ymax = ax.get_ylim()

# Load coordinates of schools in Sabah and rescale
crs={'init':'epsg:4326'}
sch = pd.read_csv('sabah/sabah_schools.csv',usecols=['lat','lon'])
sch.lat = ymin + (sch.lat - lat_min)/(lat_max - lat_min) * (ymax - ymin)
sch.lon = xmin + (sch.lon - lon_min)/(lon_max - lon_min) * (xmax - xmin)

对于 (2)，灰色背景来自于hillshade数组具有映射到灰色的地图区域之外的值。要去除灰色，我们需要取消这些值。

在这种特定情况下，我们可以利用我们知道这张地图的右上角在地图“外部”这一事实（世界上每个国家都至少有一个角落是这样的，因为没有一个国家是完美的正方形）：

top_right = hillshade[0,-1]
hillshade[hillshade == top_right] = np.nan

瞧，美丽的白色背景：

在此处输入图像描述

对于 (3)，我怀疑它需要我们以类似于我们重新缩放坐标的方式从 shapefile 重新缩放多边形。

于 2022-02-27T11:28:01.913 回答