之前写过一篇文档,是利用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))
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)
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)
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(target_ds, [1], m_layer, options=["ATTRIBUTE=RASTERVALU"])
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表示该像
元所在区域为非面要素。