【代码】多个shapefile面合并

import numpy as np
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import matplotlib.pylab as plt

# 定义一个二维数组，表示图像数据
array = np.array(( (0.1, 0.2, 0.3, 0.4),
                   (0.2, 0.3, 0.4, 0.5),
                   (0.3, 0.4, 0.5, 0.6),
                   (0.4, 0.5, 0.6, 0.7),
                   (0.5, 0.6, 0.7, 0.8) ))

# 定义两个二维数组，分别表示每个像素的纬度和经度
lat = np.array(( (10.0, 10.0, 10.0, 10.0),
                 ( 9.5,  9.5,  9.5,  9.5),
                 ( 9.0,  9.0,  9.0,  9.0),
                 ( 8.5,  8.5,  8.5,  8.5),
                 ( 8.0,  8.0,  8.0,  8.0) ))
lon = np.array(( (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5) ))

# 根据经纬度数据，计算出地理坐标范围和分辨率
xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
nrows, ncols = np.shape(array)
xres = (xmax - xmin) / float(ncols)
yres = (ymax - ymin) / float(nrows)

# 定义地理转换参数
geotransform = (xmin, xres, 0, ymax, 0, -yres)

# 创建一个新的GeoTiff文件
output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif', ncols, nrows, 1, gdal.GDT_Float32)

# 设置GeoTiff文件的地理转换参数
output_raster.SetGeoTransform(geotransform)

# 创建一个空间参考对象，设置其坐标系统为WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# 将空间参考对象的坐标系统信息导出到GeoTiff文件
output_raster.SetProjection(srs.ExportToWkt())

# 将二维数组写入GeoTiff文件
output_raster.GetRasterBand(1).WriteArray(array)

# 将所有挂起的写入操作应用到实际文件
output_raster.FlushCache()