之前写过一篇文档,是利用ArcGIS完成的本文的操作,但过程十分繁琐,其中主要的一个过程在于点矢量转栅格的过程,因此受文章 python实现使用GDAL实现矢量转栅格 的启发,利用gdal.RasterizeLayer进行矢量转栅格的操作后,实现本文的要求。
代码如下所示:

from osgeo import gdal, gdalconst
from osgeo import ogr
import gdalconst
import os
from os.path import join
import numpy as np
raster_root = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Temperature Cut'  # 原影像
shp_root = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Dian to shp 1'  # 裁剪矩形
save_root = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Dian to shp result'
raster_file = []
shp_file = []
raster_names = sorted(os.listdir(raster_root))
# 得到所需的shp文件及tif栅格图像
for raster_name in raster_names:
    raster_nameend = os.path.splitext(raster_name)[-1]
    if (raster_nameend == '.tif'):
        raster_file.append(raster_name)
shp_names = sorted(os.listdir(shp_root))
for shp_name in shp_names:
    shp_nameend = os.path.splitext(shp_name)[-1]
    if (shp_nameend == '.shp'):
        shp_file.append(shp_name)
# 理论上一个矢量对应一张栅格,因此用栅格的数量来做循环的判断
for i in range(len(raster_file)):
    raster_for = os.path.splitext(raster_file[i])[0]
    shp_for = os.path.splitext(shp_file[i])[0]
    save_name = shp_for + 'Dian.tif'
    save_path = join(save_root, save_name)
    # 得到栅格文件的最终路径
    rasterFile = join(raster_root, raster_file[i])
    shpFile = join(shp_root, shp_file[i])
    dataset = gdal.Open(rasterFile, gdalconst.GA_ReadOnly)
    raster_data = dataset.GetRasterBand(1).ReadAsArray().astype(np.float32)
    # GetGeoTransform:得到左上角横坐标、纵坐标,水平空间及垂直空间分辨率等内容
    geo_transform = dataset.GetGeoTransform()
    cols = dataset.RasterXSize  # 列数
    rows = dataset.RasterYSize  # 行数
    x_min = geo_transform[0]
    y_min = geo_transform[3]
    pixel_width = geo_transform[1]
    # 打开矢量文件
    shp = ogr.Open(shpFile, 0)
    # GetLayerByIndex:获取图层个数,shp文件一般为一层
    m_layer = shp.GetLayerByIndex(0)
    # 实现矢量转换为栅格
    target_ds = gdal.GetDriverByName('GTiff').Create(save_path, xsize=cols, ysize=rows, bands=1,
                                                     eType=gdal.GDT_Float32)
    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(dataset.GetProjection())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(0)
    band.FlushCache()
    # gdal.RasterizeLayer:根据shp或geojson矢量创建栅格
    # 跟shp字段给栅格像元赋值(本次试验中是RASTERVALU字段的值)
    gdal.RasterizeLayer(target_ds, [1], m_layer, options=["ATTRIBUTE=RASTERVALU"])
    # gdal.RasterizeLayer(target_ds, [1], m_layer) # 多边形内像元值的全是255
    del dataset
    del target_ds
    shp.Release()
    # 重新打开矢量转栅格的结果
    ShpToDian_file = gdal.Open(save_path, gdalconst.GA_ReadOnly)
    ShpToDian_data = ShpToDian_file.GetRasterBand(1).ReadAsArray().astype(np.float32)
    # 用原始栅格图像减去矢量转栅格的结果,以此来删除目标像元
    differance = raster_data - ShpToDian_data
    differance = differance.astype(np.float32)
    differance[differance == 0] = np.nan
    # 删除结果保存的文件名及路径设置
    differance_savename = 'R' + shp_for.split('H')[0] + 'T.tif'
    differance_saveroot = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Temperature Result'
    differance_savepath = join(differance_saveroot, differance_savename)
    # 输出数据
    # 设置坐标系等
    gtiff_driver = gdal.GetDriverByName("GTiff")
    out_ds = gtiff_driver.Create(differance_savepath,
                                 differance.shape[1], differance.shape[0], 1, gdal.GDT_Float32)
    # 设置输出数据坐标投影为原始影像的坐标投影
    out_ds.SetProjection(ShpToDian_file.GetProjection())
    out_ds.SetGeoTransform(ShpToDian_file.GetGeoTransform())
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(differance)
    out_band.FlushCache()
                    之前写过一篇文档,是利用ArcGIS完成的本文的操作,但过程十分繁琐,其中主要的一个过程在于点矢量转栅格的过程,因此受文章python实现使用GDAL实现矢量转栅格的启发,利用gdal.RasterizeLayer进行矢量转栅格的操作后,实现本文的要求。代码如下所示:from osgeo import gdal, gdalconstfrom osgeo import ogrimport gdalconstimport osfrom os.path import joinimport numpy
				
如图,想把背景值去掉,办法是复制栅格 工具箱>>数据管理工具>>栅格>>栅格数据集>>复制栅格 在复制栅格的窗口,在忽略背景值和NoData值处填上你的想要删除的背景值或者特定值即可。 Finish!
今天在做毕设的时候遇到栅格擦除的问题,百度了好一会儿才找到解决方案。目前在ArcMap里没有直接栅格擦除的工具,就只能自己想办法。 问题:如何在南昌栅格数据里去掉水域矢量数据? 方案:1)水域矢量水域栅格; 2)水域栅格置零; 3)将南昌栅格与水域栅格重叠部分置零; 4)南昌栅格去零。 工具:用到【面栅格】和三...
ArcMap使用指南——获取矢量空间位置对应栅格数据值简介数据步骤 平时的应用,我们在ArcMap打开了两个图层,其一个为矢量图层,一个为栅格图层,二者具有相同的空间范围,我们有时候需要知道每个要素对应栅格值。笔者第一次十分愚蠢的一个一个的使用Identify查询每个要素对应栅格值,实际上ArcMap提供了十分便捷的工具来实现这一功能。 这里用到的矢量...
3. 在“换工具”窗口,选择“矢量数据栅格数据”工具,然后击“打开”。 4. 在“矢量数据栅格数据”对话框,选择需要换的矢量数据图层,设置输出栅格数据的路径和名称。 5. 设置栅格格大小、栅格范围等参数,根据需要选择其他参数,然后击“确定”按钮。 6. 等待换过程完成,即可在指定的输出路径找到换后的栅格数据。 需要注意的是,在进行矢量数据栅格数据的过程,需要根据具体的数据类型和换要求来设置参数。例如,对于面要素,可以选择将其换为栅格数据的像值为1或0,在换后得到的栅格数据,像值为1表示该像所在区域为面要素,像值为0表示该像所在区域为非面要素。