如何用python做出一个对应shapfile文件的格网

这篇文章的目的是告诉大家怎么用python批量做一个shapefile内的规则格网。

其实arcgis已经能够完成这些问题，但为了更省事，我们可以用一些程序来解决他们。

（一）不多说，直接上代码吧！！

import os
import numpy as np
import math
import pandas as pd
from tqdm import tqdm
from shapely.geometry import Polygon, Point, MultiPolygon
root_path=os.getcwd()
shp_dir=os.path.join(root_path,'shp')
shp_paths=[os.path.join(shp_dir,i) for i in os.listdir(shp_dir)if '.shp'in i[-4:]]
result_dir=os.path.join(root_path,'grid')
#------------------------------------------------------------------------------
grid_size_list=[0.5,1]
#grid_size_list=[0.05]
def get_rec_grid(coords,grid_size_list):
    grid_dic={}
    left_lon   =math.floor(coords[0])-6
    bottom_lat =math.floor(coords[1])-6
    right_lon  =math.floor(coords[2])+6
    top_lat    =math.floor(coords[3])+6
    for i_size in grid_size_list:
        lon_grid=[i_grid for i_grid in np.arange(left_lon+i_size/2,right_lon,i_size)]
        lat_grid=[i_grid for i_grid in np.arange(bottom_lat+i_size/2,top_lat,i_size)]
        lon_lat=pd.DataFrame([[i_lon,i_lat]for i_lon in lon_grid for i_lat in lat_grid],columns=['lon','lat'])
        grid_dic.update({str(i_size):lon_lat})
    return grid_dic
def judge_point_in_shp(point_df,shp):
    
    grid_rectangle=point_df
    polygon = shapefile.Reader(shp)
    polygon = polygon.shapes() 
    shpfilePoints = [shape.points for shape in polygon]
    polygons = shpfilePoints
    point_in_list=[]
    for polygon in tqdm(polygons):
        poly = Polygon(polygon)           
        for i in range(np.shape(grid_rectangle)[0]):
            point = Point(grid_rectangle.iloc[i,0],grid_rectangle.iloc[i,1])
            if poly.contains(point):
                point_in_list.append((grid_rectangle.iloc[i,0],grid_rectangle.iloc[i,1]))
    print('over')
    return point_in_list

#------------------------------------------------------------------------------
import shapefile
import fiona
for i_shp in shp_paths:
    grid_path=os.path.join(result_dir,i_shp.split('\\')[-1].strip('.shp'))
    if os.path.exists(grid_path):
        continue
    else:
        os.mkdir(grid_path)
    shp = fiona.open(i_shp)
    coords = shp.bounds
    shp.close()
    rec_grid_dic=get_rec_grid(coords,grid_size_list)
    for i_key in rec_grid_dic.keys():
        point_df=rec_grid_dic[i_key]
        result=judge_point_in_shp(point_df,i_shp)#+judge_point_in_buffershp(point_df,i_shp,float(i_key))
        result=pd.DataFrame(result,columns=['lon','lat'])
        result.drop_duplicates()
        result.to_csv(os.path.join(grid_path,i_shp.split('\\')[-1].strip('.shp')+'_'+i_key+'.csv'),index=None)
        

（二）代码使用中需要注意的点：

①如果程序报错，可以先检查一下自己的shapefile文件的属性表中是不是有中文。

②文件夹的相对位置

 ③这种方法由于涉及到判断点在shp的polygon内与否，所以当点非常多的时候不太适合用这种方法，比如提取中国的0.1°格网的话可能会耗时比较多，这种情况可以用arcgis提取了

④代码中有判断格网文件夹是否存在的功能，如果已存在但没有提取成功，再次运行需要删掉相对应的文件夹。 

⑤格网有偏移的话，需要调整get_rec_grid函数。

⑥在实际的应用中，离polygon非常近的点有的时候也是需要的，可以用arcgis加做个buffer

⑦有啥问题自己先尝试一下，看看文件夹啥的，改改路径或者名称等，程序很简单。。。