快捷導(dǎo)航

利用Python柵格化地圖(以成都市為例,含代碼)

更新時間：2024年03月14日 09:00:27 作者：數(shù)據(jù)的旅途

這篇文章主要給大家介紹了關(guān)于利用Python柵格化地圖的相關(guān)資料,
Python中可以使用多種庫來進(jìn)行柵格化地圖的操作,其中比較常用的有g(shù)eopandas、rasterio等,文中通過代碼介紹的非常詳細(xì),需要的朋友可以參考下

python代碼實(shí)現(xiàn)

讀取成都市邊界的圖層文件(.shp)，并可視化

import geopandas as gpd

cd_shape = gpd.read_file('Chengdu/Chengdu.shp')
cd_shape.plot(edgecolor='k',facecolor='none')

下面這個柵格地圖的類是我自己寫的，類的參數(shù)主要有

rasterDataPath：圖層文件的路徑
length：柵格單元的長度，單位為 m

class RasterData:

    def __init__(self, raster_data_path: str, length: float):
        self.raster = gpd.read_file(raster_data_path)
        self.length = length/1000 * 0.009
        self.polygons = []
        self.grid_ids = []
        self.x_min, self.y_min, self.x_max, self.y_max = self.raster.total_bounds
        self.rows, self.cols = self.grid_shape()

    def grid_shape(self) -> tuple:
        rows = int(math.ceil((self.y_max - self.y_min) / float(self.length)))
        cols = int(math.ceil((self.x_max - self.x_min) / float(self.length)))
        return rows, cols

    def grid_num(self) -> int:
        return self.rows * self.cols

    def grid_map(self) -> gpd.GeoDataFrame:
        points_list = []
        for row in range(self.rows):
            for col in range(self.cols):
                center_point_x = self.x_min + self.length / 2 + col * self.length
                center_point_y = self.y_min + self.length / 2 + row * self.length
                
                points = [Point(center_point_x + dx * self.length / 2,
                                center_point_y + dy * self.length / 2)
                          for dx, dy in [(-1, 1), (1, 1), (1, -1), (-1, -1)]]
                points_list.append(points)

        polygons = [Polygon(points) for points in points_list]
        grid_ids = list(range(len(polygons)))
        
        grid = gpd.GeoDataFrame({'geometry': polygons, 'grid_id': grid_ids}, crs=self.raster.crs)
        return grid
        
    # 計算柵格與區(qū)域的交集
    def grid_intersection(self, region: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
        grid = self.grid_map()
        intersection_data = gpd.overlay(grid, region, how='intersection')
        return intersection_data

實(shí)例化對象并調(diào)用grid_map方法

grid = RasterData('Chengdu/Chengdu.shp', 2000) # 實(shí)例化對象
grid_data = grid.grid_map() # 調(diào)用grid_map方法進(jìn)行柵格化
# 將兩個圖層繪制在一起
fig, ax = plt.subplots(figsize=(10, 10))
# 加粗繪圖的線寬
cd_shape.plot(ax=ax, edgecolor='k', linewidth=1, facecolor='none')
grid_data.plot(ax=ax, edgecolor='k', linewidth=0.5, facecolor='none')

得到柵格模型，但此時的柵格是根據(jù)成都市邊界的最大范圍進(jìn)行劃分的，很多時候我們需要的是地理邊界內(nèi)部的柵格，因此需要調(diào)用grid_intersection方法

# 將兩個圖層繪制在一起
fig, ax = plt.subplots(figsize=(10, 10))
# 加粗繪圖的線寬
intersection_data = grid.grid_intersection(cd_shape)
cd_shape.plot(ax=ax, edgecolor='k', linewidth=1, facecolor='none')
intersection_data.plot(ax=ax, edgecolor='k', linewidth=0.5, facecolor='none')