【Arcpy】提取遥感影像像元值

燕南路GISer

已于 2022-08-12 22:12:21 修改

阅读量2.6k

点赞数 1

分类专栏： Arcpy 文章标签： python 开发语言

于 2022-08-12 17:07:48 首次发布

本文链接：https://blog.youkuaiyun.com/xza13155/article/details/126278899

版权

Arcpy 专栏收录该内容

19 篇文章

订阅专栏

基于开源库执行

GDAL 是读写大量的栅格空间数据格式的广泛应用的开源库。该库起源于 1998年,已经大幅进化。它支持他自己的数据模型和应用程序接口(API)。从最初的单一发展的起源,GDAL已发展成为一个分布式的项目,开发人员的数量相对比较大。
目前GDAL是OSGEO的子项目，代码进行了重新组织。在 RFC 17号文件中，实现了Python的新的名称空间osgeo，并将gdal与ogr都包含在这个名称空间之下。

下面给出主要学习的链接:
官方文档
 python与开源GIS
OSGeo中国中心

下面附上基于gdal计算像元值的代码。

import os
import math
import numpy as np
from tqdm import tqdm
from openpyxl import Workbook
from osgeo import gdal, gdalconst
#读取python文件所在的位置
base_dir = os.path.dirname(os.path.abspath(__file__))


#读取图像数据
def getData(path):
    #读取一个GeoTiff文件的信息，GA_ReadOnly只读形式
    dataset = gdal.Open(path, gdalconst.GA_ReadOnly)
    #栅格数据的宽度(X方向上的像素个数)
    width = dataset.RasterXSize
    #栅格数据的高度(Y方向上的像素个数)
    height = dataset.RasterYSize
    
    '''
六个参数分别是：
    geos[0]  top left x 左上角x坐标
    geos[1]  w-e pixel resolution 东西方向像素分辨率
    geos[2]  rotation, 0 if image is "north up" 旋转角度，正北向上时为0
    geos[3]  top left y 左上角y坐标
    geos[4]  rotation, 0 if image is "north up" 旋转角度，正北向上时为0
    geos[5]  n-s pixel resolution 南北向像素分辨率，为负值
    x/y为图像的x/y坐标，geox/geoy为对应的投影坐标
    
    (294256.45800898736, 500.0, 0.0, 4314092.674585257, 0.0, -500.0)
    可知，200201_month_p的像元大小为500
'''
    #栅格数据仿射变换的六参数
    im_geotrans = dataset.GetGeoTransform()
    #图像的投影信息
    im_proj = dataset.GetProjection()
    
    '''
    ReadAsArray(xoff=0, yoff=0, xsize=None, ysize=None)
    xoff,yoff ：指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位）。
    xsize,ysize ： 指定要读取部分图像的矩形的长和宽（以像元为单位）。
    
    dataset.ReadAsArray(0, 0, width, height)
    就是读取全图
    '''
    
    #读取图像数据(以numpy数组的形式)
    data = np.array(dataset.ReadAsArray(0, 0, width, height), dtype=float)
    #返回图像数据，六参数和投影
    return data, im_geotrans, im_proj

#换算在array中的下标
def get_r_c(lat, lon, im_geotrans):
    row = math.floor((im_geotrans[3] - lat) / np.abs(im_geotrans[5]))
    col = math.floor((lon - im_geotrans[0]) / im_geotrans[1])
    return row, col
#写入excel
def get_ref_excel(data, im_geotrans):
    wbook = Workbook()
    book = wbook.create_sheet(0)
    for i in range(data.shape[0]):
        lat = im_geotrans[3] + i*im_geotrans[5]
        col = 2
        for j in range(data.shape[1]):
            lon = im_geotrans[0] + j*im_geotrans[1]
            
            book.cell(row=i+2, column=1).value = lat
            book.cell(row=1, column=col).value = lon
            book.cell(row=i+2, column=col).value = data[i,j]
            col += 1
    wbook.save(os.path.join(base_dir, '201812_sm.xlsx'))
    
    
#输出modis数据
def writeRSData(filename, im_data, im_geotrans, im_proj): #将数据写入Tif文件
    path, name = os.path.split(filename)
    if not os.path.exists(path):
        os.makedirs(path)
        
    if 'int8' in im_data.dtype.name:  # 判断栅格数据的数据类型
            datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
            
    if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape  # 多维或1.2维
        
    driver = gdal.GetDriverByName("GTiff")            #数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
    dataset.SetGeoTransform(im_geotrans)              #写入仿射变换参数
    dataset.SetProjection(im_proj)                    #写入投影
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  #写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])
    del dataset
    
#读取图像
#ref_file, im_geotrans_ref, im_proj_ref分别为图像数据，六参数和投影
ref_file, im_geotrans_ref, im_proj_ref = getData(r'data\\200201_month_p.tif')
extract_file, im_geotrans_c, im_proj_c = getData(r'data\\200201_1month_NDVI.tif')

#创建一个excel表
#sheet工作表的集合
wbook = Workbook()
# row 行 或 column 列的一维集合
book = wbook.create_sheet(0)

#列表化图像的仿射变换参数
im_geotrans_ref = list(im_geotrans_ref)
im_geotrans_c = list(im_geotrans_c)
print(im_geotrans_ref)

#get_ref_excel(ref_file, im_geotrans_ref) 输出modis数据，不需要的话就注释掉

#tqdm是 Python 进度条库，可以在 Python长循环中添加一个进度提示信息。
#ref_file.shape[0]代表行数，ref_file.shape[1]代表列数
#遍历行（经度）
for i in tqdm(range(ref_file.shape[0]), desc='提取数据'):
    #读取经度(左上角初始经度，每次偏移一个像素点)
    lat = im_geotrans_ref[3] + i*im_geotrans_ref[5]
    col = 2
    #遍历列（维度）
    for j in range(ref_file.shape[1]):
        #读取经度(左上角初始维度，每次偏移一个像素点)
        lon = im_geotrans_ref[0] + j*im_geotrans_ref[1]
        #计算出array中对应的值的位置
        row_c = math.floor(np.abs(im_geotrans_c[3] - lat) / np.abs(im_geotrans_c[5]))
        col_c = math.floor(np.abs(lon - im_geotrans_c[0]) / np.abs(im_geotrans_c[1]))
        #填写维度，第一列的那些数据
        book.cell(row=i+2, column=1).value = lat
        #填写经度，第一行的数据
        book.cell(row=1, column=col).value = lon
        #填写经纬度对应的值
        book.cell(row=i+2, column=col).value = extract_file[row_c, col_c]
        col += 1
print(base_dir)
wbook.save(os.path.join(base_dir, '200201_ndvi_1month.xlsx'))

导出为图片

# coding=GBK
from osgeo import gdal
import pandas
import numpy as np

csv_data = pandas.read_csv('200201_ndvi_month.csv')#CSV所在位置
np_data = np.array(csv_data)
print(np_data)
#np_data[np_data == -9999] = np.nan
b4dataset = gdal.Open('200201_month_p.tif')#投影信息——将现有投影信息赋给将要生成的tif
im_geotrans = b4dataset.GetGeoTransform()  # 获取仿射矩阵信息
im_proj = b4dataset.GetProjection()  # 获取投影信息
datatype = gdal.GDT_Float64
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create('200201_ndvi_month.tif', np_data.shape[1], np_data.shape[0], 1, datatype)
#tif影像存放位置
dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
dataset.SetProjection(im_proj)  # 写入投影
dataset.GetRasterBand(1).WriteArray(np_data)

基于arcpy

将要提取的坐标转为点数据，调用arcpy的值提取到点工具执行，然后将结果导出为表格

# coding=GBK
#读取坐标
lines = open('point_inform.csv','r')
i=1
for line in lines:
    if(i==1):
        lats=line.strip('\n').strip(',').split(',')
    else:
        lons=line.strip('\n').split(',')
    i=i+1
#将坐标输出为csv文件
fo = open('point.csv','w')
fo.writelines('ID,lon,lat\n')
fo.close()
i=1
fo = open('point.csv','a')
for lat in lats:
    for lon in lons:
        line=str(i)+','+lon+','+lat
        fo.writelines(line+'\n')
        i=i+1
fo.close()

# coding=GBK

import arcpy
from arcpy import env
arcpy.CheckOutExtension('Spatial')
data_path="pythonProject\\"
image_path="pythonProject\data\\"
arcpy.env.workspace = "pythonProject\\数据库.gdb"
#可以覆盖已有要素
env.overwriteOutput = True
'''
#转点要素
arcpy.MakeXYEventLayer_management(data_path+"point.csv", "lon", "lat", "point_t","WGS 1984 UTM Zone 50N.prj")
arcpy.FeatureToPoint_management ("point_t","point")
#从NDVI提取值
arcpy.sa .ExtractValuesToPoints ("point",image_path+'200201_1month_NDVI.tif',"point_value")
'''
#获取经纬度信息
lines = open(data_path+'point_inform.csv','r')
i=1
for line in lines:
    #维度
    if(i==1):
        lats=line.strip('\n').strip(',').split(',')
        lats_len=len(lats)
    #经度
    else:
        lons=line.strip('\n').strip(',').split(',')
        lons_len = len(lons)
    i=i+1
lines.close()
#写入其余行
cursor = arcpy.SearchCursor("point_value")
fo = open('200201_ndvi_month.csv','w')
n=0
i=0
line=''
for row in cursor:
    if(line==''):
        line= str(row.getValue("RASTERVALU"))
    else:
        line=line +','+ str(row.getValue("RASTERVALU"))
    n=n+1
    if(n==lons_len):
        fo.writelines(line+'\n')
        n=0
        i=i+1
        if(i==lats_len):
            break
        line=''
fo.close()
'''
#写入表格第一行
fo = open('200201_ndvi_month.csv','w')
line=' '
for lon in lons:
    line=line+','+lon
fo.writelines(line+'\n')
fo.close()
#写入其余行
cursor = arcpy.SearchCursor("point_value")
fo = open('200201_ndvi_month.csv','a')
n=1
i=0
line=lats[i]
for row in cursor:
    line=line +','+ str(row.getValue("RASTERVALU"))
    n=n+1
    if(n==lons_len):
        fo.writelines(line+'\n')
        n=1
        i=i+1
        if(i==lats_len):
            break
        line=lats[i]
fo.close()
'''