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的值提取到点工具执行，然后将结果导出为表格
# 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()
'''