Python 使用 GDAL 和 geopandas 处理全国建成区数据

使用Python中的GDAL和geopandas库对全国建成区数据（tif格式）进行处理，将每个栅格点赋值为行政区划编号和原始值的组合，实现建成区数据与行政区划数据的关联。\n\n首先，需要安装GDAL和geopandas库。可以使用以下命令进行安装：\n\n\npip install GDAL\npip install geopandas\n\n\n接下来，可以使用以下代码来实现对tif文件的处理：\n\npython\nimport gdal\nimport geopandas as gpd\nimport numpy as np\n\n# 读取tif文件\ntif_path = 'path/to/your/tif/file.tif'\nds = gdal.Open(tif_path)\nband = ds.GetRasterBand(1)\narray = band.ReadAsArray()\n\n# 读取shp文件\nshp_path = 'path/to/your/shp/file.shp'\ndf = gpd.read_file(shp_path)\n\n# 对栅格重新赋值\nfor index, row in df.iterrows():\n geom = row.geometry\n mask = geom.contains(gpd.points_from_xy(ds.GetGeoTransform()[0] + ds.GetGeoTransform()[1] * np.arange(array.shape[1]),\n ds.GetGeoTransform()[3] + ds.GetGeoTransform()[5] * np.arange(array.shape[0])))\n array[mask] = index * 10 + array[mask]\n\n# 保存处理后的tif文件\nnew_tif_path = 'path/to/save/processed/tif/file.tif'\ndriver = gdal.GetDriverByName('GTiff')\nnew_ds = driver.Create(new_tif_path, ds.RasterXSize, ds.RasterYSize, 1, band.DataType)\nnew_ds.SetGeoTransform(ds.GetGeoTransform())\nnew_ds.SetProjection(ds.GetProjection())\nnew_ds.GetRasterBand(1).WriteArray(array)\nnew_ds.FlushCache()\nnew_ds = None\n\nprint("处理完成")\n\n\n上述代码中，首先使用gdal库的Open函数读取tif文件，并读取第一个波段的数据。然后使用geopandas库的read_file函数读取shp文件。\n\n接下来，使用iterrows函数遍历shp文件的每一行数据，获取每个行政区划的几何对象。然后使用contains函数判断每个栅格点是否在当前行政区划内，得到一个布尔掩码。最后，根据布尔掩码将对应栅格点的值重新赋值为行政区划的编号和原始值的组合。\n\n最后，使用gdal库创建一个新的tif文件，并将处理后的栅格数据写入其中。\n\n注意，上述代码中的路径需要根据实际情况进行修改。