def plot_map(self, data, city_info, source_city): try: fig, ax = plt.subplots(figsize=(12, 10)) basemap = Basemap(llcrnrlon=city_info["llcrnrlon"], llcrnrlat=city_info["llcrnrlat"], urcrnrlon=city_info["urcrnrlon"], urcrnrlat=city_info["urcrnrlat"]) city_name = city_info["city"] # 读取 YRD 地图轮廓 basemap.readshapefile(r'..\data\shp_citys\{}'.format(city_name), city_name) # 提取所有有效多边形 polygons = [] shapes = getattr(basemap, city_name) for shape in shapes: poly = sgeom.Polygon(shape) if poly.is_valid: polygons.append(poly) if not polygons: raise ValueError("未能从 shapefile 中提取有效多边形") # 合并多边形 merged_poly = unary_union(polygons) # 根据几何类型构造 clip path if merged_poly.geom_type == 'Polygon': x, y = merged_poly.exterior.xy verts = np.array(list(zip(x, y))) clip_path = Path(verts) elif merged_poly.geom_type == 'MultiPolygon': # 对于 MultiPolygon,需要创建复合路径 paths = [] for poly in merged_poly.geoms: x, y = poly.exterior.xy verts = np.array(list(zip(x, y))) # 构造路径码:第一个点为 MOVETO,后续为 LINETO codes = [Path.MOVETO] + [Path.LINETO] * (len(verts) - 1) paths.append(Path(verts, codes)) clip_path = Path.make_compound_path(*paths) else: raise ValueError("不支持的几何类型用于裁剪") # 绘制底图数据 cm = cmaps.precip3_16lev im = ax.imshow(data, interpolation='None', extent=(70, 140, 60, 0), cmap=cm, norm=LogNorm(vmin=1, vmax=10000)) # 创建 patch 并设置 clip path patch = PathPatch(clip_path, transform=ax.transData) im.set_clip_path(patch) # 添加颜色条 divider = make_axes_locatable(ax) cax = divider.append_axes('right', size='5%', pad=0.05) plt.colorbar(im, cax=cax) output_file = f"./distribution_map/2019nox/2019nox_{source_city}_to_{city_info['city_eng_name']}.png" plt.savefig(output_file, dpi=1000) print(f"目标城市 {source_city} 对影响城市 {city_info['city_eng_name']} 的图绘制完成") except Exception as e: print(f"Error plotting map: {e}")解释一下这段代码,并详细展开介绍其中涉及的函数和用法
时间: 2025-07-07 22:47:28 浏览: 11
<think>嗯,我现在要仔细看一下用户提供的这段Python代码,并尝试理解它的功能。这段代码是一个名为`plot_map`的方法,属于某个类的一部分,因为它使用了`self`作为第一个参数。代码的主要目的是绘制地图,并保存为图片。让我一步步拆解分析。
首先,用户希望我解释这段代码,并详细展开其中涉及的函数和用法。我需要确保自己理解每一行代码的作用,以及所使用的库和函数。可能需要回忆一下Basemap库的相关知识,以及matplotlib的一些用法。
首先,函数定义部分:
```python
def plot_map(self, data, city_info, source_city):
```
这是一个类的方法,接受三个参数:data(数据)、city_info(城市信息字典)、source_city(源城市名称)。接下来是try块,说明代码可能抛出异常,并在except中处理。
在函数内部,首先创建了fig和ax对象:
```python
fig, ax = plt.subplots(figsize=(12, 10))
```
这里使用了matplotlib的subplots方法创建了一个图形和一个坐标轴对象,设置图形大小为12x10英寸。
接下来是Basemap的初始化:
```python
basemap = Basemap(llcrnrlon=city_info["llcrnrlon"],
llcrnrlat=city_info["llcrnrlat"],
urcrnrlon=city_info["urcrnrnrlon"],
urcrnrlat=city_info["urcrnrlat"])
```
Basemap是Matplotlib的一个工具库,用于绘制地图。这里用city_info中的经纬度设定了地图的边界(左下角和右上角)。需要注意的是,可能用户输入时是否有拼写错误,比如urcrnrlon是否正确,因为代码中urcrnrlon在city_info中是否对应正确的键名?比如是否应该是urcrnrlon而不是urcrnrlon?需要检查是否有拼写错误,但用户提供的代码可能存在typo,比如urcrnrlon是否正确。
然后读取shapefile:
```python
basemap.readshapefile(r'..\data\shp_citys\{}'.format(city_name), city_name)
```
这里使用了Basemap的readshapefile方法,路径是相对路径,指向存储城市形状文件的目录。city_name来自city_info中的"city"键。读取后,Basemap会将形状信息存储在属性中,可以通过`basemap.城市名`来获取,比如shapes和info。
接下来提取有效多边形:
```python
polygons = []
shapes = getattr(basemap, city_name)
for shape in shapes:
poly = sgeom.Polygon(shape)
if poly.is_valid:
polygons.append(poly)
```
这里使用getattr获取basemap对象的city_name属性,即读取的形状数据。然后遍历每个形状,将其转换为Shapely的Polygon对象,并检查是否有效(is_valid)。有效多边形被添加到列表中。这一步可能是为了过滤掉无效的地理多边形,确保后续处理正确。
合并多边形:
```python
merged_poly = unary_union(polygons)
```
这里使用了Shapely的unary_union方法,将所有有效多边形合并成一个几何对象。如果多边形有重叠,这个方法会将它们合并为一个整体,减少后续处理的复杂度。
构造裁剪路径:
根据合并后的几何类型(Polygon或MultiPolygon),生成对应的裁剪路径。如果是单个多边形,直接取外部的坐标点,构造Path对象。如果是多个多边形(MultiPolygon),则需要为每个子多边形创建路径,并使用复合路径(compound path)。
绘制数据:
```python
im = ax.imshow(data, ...)
```
使用imshow方法将数据可视化为图像,设置插值为None,颜色映射为precip3_16lev,使用对数归一化。extent参数设定了图像的经纬度范围,这里需要确认extent的值是否正确,比如(70, 140, 60, 0)可能对应经度70到140,纬度0到60,这可能覆盖中国大部分区域?需要结合city_info中的经纬度边界来看是否一致。
设置裁剪路径:
```python
patch = PathPatch(clip_path, transform=ax.transData)
im.set_clip_path(patch)
```
这里将之前生成的clip_path应用到图像上,使得只有地图轮廓内的部分显示,外部被裁剪掉。PathPatch将路径转换为Patch对象,transform指定坐标转换,确保路径与数据坐标对齐。
添加颜色条:
使用make_axes_locatable来在右侧添加一个颜色条,设置大小和间隔,然后保存图片到指定路径。
异常处理部分捕获所有异常并打印错误信息,这有助于调试但可能过于宽泛,最好能具体一些。
现在需要总结每个步骤的关键点,并解释涉及的函数和库,比如Basemap、Shapely、matplotlib的Path和PathPatch,以及LogNorm的作用。同时指出潜在问题,如路径拼写错误、坐标范围是否正确、形状文件是否存在等。
可能需要提醒用户关于路径中的反斜杠在Python字符串中需要使用原始字符串或转义,比如r'..\data...'是正确的,或者应使用正斜杠以避免转义问题。另外,读取shapefile时的路径是否正确,是否所有城市都有对应的shapefile文件,这些都是潜在的问题点。
另外,合并多边形时如果多边形无效可能引发错误,代码中已经检查了is_valid,但如果polygons为空,会抛出ValueError。用户需要确保shapefile中的多边形有效,并且至少有一个有效。
在裁剪路径构造部分,MultiPolygon的处理是否正确,是否所有子多边形都被正确包含到复合路径中,这可能影响最终的裁剪效果。
最后,保存图片的路径是否正确,是否有写入权限,以及dpi设置是否过高导致文件过大,这些也是需要考虑的。</think>以下是 `plot_map` 方法的逐层解析及关键函数说明:
---
### **代码整体功能**
这是一个基于 `Basemap` 和 `Matplotlib` 的地理数据可视化方法,核心步骤包括:
1. 读取地理边界 Shapefile
2. 构建地图裁剪路径
3. 绘制热力图并限定显示范围
4. 保存带裁剪效果的专题地图
---
### **关键步骤详解**
#### **1. 地图初始化**
```python
basemap = Basemap(llcrnrlon=city_info["llcrnrlon"],
llcrnrlat=city_info["llcrnrlat"],
urcrnrlon=city_info["urcrnrlon"],
urcrnrlat=city_info["urcrnrlat"])
```
- **作用**:创建 Basemap 实例,定义地图边界
- **参数说明**:
- `llcrnrlon/lat`:左下角经纬度 (Lower Left Corner)
- `urcrnrlon/lat`:右上角经纬度 (Upper Right Corner)
- **注意**:需确保 `city_info` 字典包含这四个键值
---
#### **2. 读取 Shapefile**
```python
basemap.readshapefile(r'..\data\shp_citys\{}'.format(city_name), city_name)
```
- **作用**:加载城市边界矢量数据
- **参数**:
- 路径格式:`..\data\shp_citys\城市名`(需包含 `.shp` 文件)
- 名称标签:`city_name` 会作为属性名存储形状数据
- **后续访问**:通过 `basemap.城市名` 获取坐标点集合
---
#### **3. 构建有效多边形**
```python
poly = sgeom.Polygon(shape)
if poly.is_valid:
polygons.append(poly)
```
- **库依赖**:`shapely.geometry`(需安装)
- **关键方法**:
- `Polygon()`:将坐标点转为几何多边形
- `is_valid`:检查几何有效性(如无自相交)
---
#### **4. 合并多边形**
```python
merged_poly = unary_union(polygons)
```
- **作用**:合并多个多边形为一个几何体(处理飞地等复杂边界)
- **原理**:使用 `shapely.ops.unary_union` 进行空间合并操作
---
#### **5. 创建裁剪路径**
```python
clip_path = Path(verts) # 单多边形
# 或多边形复合路径
clip_path = Path.make_compound_path(*paths)
```
- **相关类**:
- `matplotlib.path.Path`:定义矢量路径
- `PathPatch`:将路径转换为可绘制的图形元素
- **MOVETO/LINETO**:路径绘制指令(类似 SVG)
---
#### **6. 数据可视化与裁剪**
```python
im = ax.imshow(data, ..., norm=LogNorm(vmin=1, vmax=10000))
im.set_clip_path(patch)
```
- **关键参数**:
- `extent=(70, 140, 60, 0)`:图像映射到地图的经纬度范围(需与 Basemap 范围匹配)
- `LogNorm`:对数颜色标准化(需 `matplotlib.colors`)
- **裁剪原理**:通过 `set_clip_path` 限制图像仅在多边形区域内显示
---
#### **7. 颜色条与输出**
```python
plt.colorbar(im, cax=cax)
plt.savefig(output_file, dpi=1000)
```
- **`make_axes_locatable`**:动态创建与主图对齐的附属坐标轴
- **dpi=1000**:高分辨率输出(注意文件大小)
---
### **潜在问题与优化建议**
1. **路径风险**:
- 使用原始字符串 `r''` 防止转义错误
- 建议改用 `os.path.join` 处理文件路径
2. **坐标系一致性**:
- 确保 `imshow(extent=...)` 与 Basemap 边界对齐
- 推荐通过 `basemap.transform` 转换坐标
3. **异常处理**:
- 当前 `except Exception` 过于宽泛,建议细化异常类型(如 `KeyError` 检查字典键)
4. **性能优化**:
- 高 dpi 输出可能导致内存问题,可适当降低或分块处理
---
### **典型应用场景**
- 污染物扩散可视化(如代码中的 NOx 分布)
- 区域经济指标空间分布
- 气象数据(降水量、温度)区域渲染
是否需要进一步解释某个具体技术细节?
阅读全文
相关推荐


















