【GDAL学习】用GDAL读取栅格数据

1.根据坐标读取遥感影像的单个像素值

# week 4: get pixel values at a set of coordinates by reading in one pixel at a time
import os, sys, time
from osgeo import gdal
from gdalconst import *# start timing
startTime = time.time()# coordinates to get pixel values for
xValues = [447520.0, 432524.0, 451503.0]
yValues = [4631976.0, 4608827.0, 4648114.0]# set directory
os.chdir('E:/data/GDAL/ospy_data4')
# register all of the drivers
gdal.AllRegister()# open the image
ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:print('cannot open this file.')sys.exit(1)# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount# get georeference info
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]# loop through the coordinates
for i in range(3):# get x,yx = xValues[i]y = yValues[i]# compute pixel offsetxOffset = int((x - xOrigin) / pixelWidth)yOffset = int((y - yOrigin) / pixelHeight)# create a string to print outs = str(x) + ' ' + str(y) + ' ' + str(xOffset) + ' ' + str(yOffset) + ' '# loop through the bandsfor j in range(bands):band = ds.GetRasterBand(j + 1)# read data and add the value to the stringdata = band.ReadAsArray(xOffset, yOffset, 1, 1)value = data[0, 0]s = s + str(value) + ' '# print out the data stringprint(s)
# figure out how long the script took to run
endTime = time.time()
print('The script took ' + str(endTime - startTime) + ' seconds.')

2.根据坐标读取遥感影像的单个像素值，通过获取图像所有像素值获得

# week 4:script to get pixel values at a set of coordinates by reading in entire bands
from osgeo import gdal
from gdalconst import *
import os, sys, time# start timing
startTime = time.time()# coordinates to get pixel values for
xValues = [447520.0, 432524.0, 451503.0]
yValues = [4631976.0, 4608827.0, 4648114.0]# set directory
os.chdir('E:/data/GDAL/ospy_data4')# register all of the drivers
gdal.AllRegister()# open the image
ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:print('Cannot open aster.img')sys.exit(1)# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount# get georeference info
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]# create a list to store band data in
bandList = []# read in bands and store all the data in bandList
for i in range(bands):band = ds.GetRasterBand(i + 1)data = band.ReadAsArray(0, 0, cols, rows)bandList.append(data)# loop through the coordinates
for i in range(3):x = xValues[i]y = yValues[i]# compute pixel offsetxOffset = int((x - xOrigin) / pixelWidth)yOffset = int((y - yOrigin) / pixelHeight)# create a string to print outs = str(x) + ' ' + str(y) + ' ' + str(xOffset) + ' ' + str(yOffset) + ' '# loop through the bands and get the pixel valuefor j in range(bands):data = bandList[j]value = data[xOffset-1, yOffset-1]s = s + str(value)+' '# print out the data stringprint(s)# figure out how long the script took to run
endTime = time.time()
print('The script took ' + str(endTime - startTime) + ' seconds.')

3.Assignment 4a：Read pixel values from an image

Print out the pixel values for all three bands of aster.img at the points contained in sites.shp •
Use any method of reading the raster data that you want, but I would suggest one pixel at a time (fastest in this case since we don't need much data)

1）自己写的代码，现获取所有sities的坐标后再获取全部的pixel values

# Assignment 4a: Read pixel values from an image
# write by myself
from osgeo import gdal
from gdalconst import *
import ogr, os, sys# set the directory
os.chdir('E:/data/GDAL/ospy_data4')# set the driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# open the file
siteDS = driver.Open('sites.shp', 0)
if siteDS is None:print('Cannot open sites.shp')sys.exit(1)# get the layer
siteLy = siteDS.GetLayer()
# get the feature
siteFe = siteLy.GetNextFeature()# create lists to store the x,y coordination of the feature
sitePosX = []
sitePosY = []# cycle the feature to get the x,y coordination
while siteFe:geom = siteFe.GetGeometryRef()fex = geom.GetX()fey = geom.GetY()sitePosX.append(fex)sitePosY.append(fey)siteFe.Destroy()siteFe = siteLy.GetNextFeature()# create a txt file to store result
txt = open('output.txt', 'w')# register all of the drivers
gdal.AllRegister()
# open the image
asterDS = gdal.Open('aster.img', GA_ReadOnly)
if asterDS is None:print('Cannot open aster.img')sys.exit(1)# get the image size
rows = asterDS.RasterYSize
cols = asterDS.RasterXSize
bands = asterDS.RasterCount# get the georeference info
transform = asterDS.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]# get the pixel value of the coordination
for i in range(len(sitePosX)):x = sitePosX[i]y = sitePosY[i]# get the offsetxOffset = int((x - xOrigin) / pixelWidth)yOffset = int((y - yOrigin) / pixelHeight)s = str(i + 1) + ' ' + str(x) + ' ' + str(y) + ' ' + str(xOffset) + ' ' + str(yOffset) + ' '# get all bands of the pixel valuefor j in range(bands):band = asterDS.GetRasterBand(j + 1)data = band.ReadAsArray(xOffset, yOffset, 1, 1)value = data[0, 0]s = s + str(value) + ' 's = s + '\n'# write to a text filetxt.write(s)print(s)# destroy datasource
siteDS.Destroy()

2）参考答案给的代码，每获取一个site坐标就获取该坐标的pixel value

import os, sys
try:from osgeo import ogr, gdalfrom osgeo.gdalconst import *os.chdir('E:/data/GDAL/ospy_data4')
except ImportError:import ogr, gdalfrom gdalconst import *os.chdir('E:/data/GDAL/ospy_data4')# open the shapefile and get the layer
driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open('sites.shp')
if shp is None:print('Could not open sites.shp')sys.exit(1)
shpLayer = shp.GetLayer()# register all of the GDAL drivers
gdal.AllRegister()# open the image
img = gdal.Open('aster.img', GA_ReadOnly)
if img is None:print('Could not open aster.img')sys.exit(1)# get image size
rows = img.RasterYSize
cols = img.RasterXSize
bands = img.RasterCount# get georeference info
transform = img.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]# loop through the features in the shapefile
feat = shpLayer.GetNextFeature()
while feat:geom = feat.GetGeometryRef()x = geom.GetX()y = geom.GetY()# compute pixel offsetxOffset = int((x - xOrigin) / pixelWidth)yOffset = int((y - yOrigin) / pixelHeight)# create a string to print outs = feat.GetFieldAsString('ID') + ' '# loop through the bandsfor j in range(bands):band = img.GetRasterBand(j + 1)  # 1-based index# read data and add the value to the stringdata = band.ReadAsArray(xOffset, yOffset, 1, 1)value = data[0, 0]s = s + str(value) + ' '# print out the data stringprint(s)# get the next featurefeat.Destroy()feat = shpLayer.GetNextFeature()# close the shapefile
shp.Destroy()

4.为了更有效率地读取栅格，一种更好的方法是分块读取数据。数据文件中给出了一个utils.py文件中的GetBlockSize()函数可以实现获取图像中的xBlockSize和yBlockSize，但是在我的电脑里不能运行。于是尝试手动定义xBlockSize和yBlockSize，更改了几次数值发现程序运行结果都是一样的，因此我认为手动赋值给xBlockSize和yBlockSize也是可行的。下面代码是循环分块读取栅格数据的样例，具体原理可参见PPT

from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys
import utilsos.chdir('E:/data/GDAL/ospy_data4')
gdal.AllRegister()ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:print('Cannot open aster.img')sys.exit(1)rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCountband = ds.GetRasterBand(1)
# 由于utils.py文件中找不到_gdal模块，因此以下三行代码不能使用
# xblocksize和yblocksize改为手动设置
# 更改了几个值发现程序运行最终结果不变
# blockSizes = utils.GetBlockSize(band)
# xBlockSize = blockSizes[0]
# yBlockSize = blockSizes[1]xBlockSize = 20
yBlockSize = 20count = 0for i in range(0, rows, yBlockSize):if (i + yBlockSize) < rows:numRows = yBlockSizeelse:numRows = rows - ifor j in range(0, cols, xBlockSize):if (j + xBlockSize) < cols:numCols = xBlockSizeelse:numCols = cols - jdata = band.ReadAsArray(j, i, numCols, numRows).astype(np.float)mask = np.greater(data, 0)count = count + np.sum(np.sum(mask))print('Number of non-zero pixels: ', count)

5. Assignment 4b：

Write a script to calculate the average pixel value for the first band in aster.img
Read in the data one block at a time
Do the calculation two ways ：
- Including zeros in the calculation
- Ignoring zeros in the calculation

1）自己写的代码

# Assignment 4b:Write a script to calculate the average pixel value for the first band in aster.img
# write by myself
# import models
from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys# set the directory
os.chdir('E:/data/GDAL/ospy_data4')# register all drivers
gdal.AllRegister()# open the file
ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:print('Cannot open aster.img')sys.exit(1)# get the raster size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount# get the first band and set blocksize
band = ds.GetRasterBand(1)
xBlockSize = 10
yBlockSize = 10sumIncludeZero = 0.0
zeroCount = 0# cycle row blocksize
for i in range(0, rows, yBlockSize):# set yblocksizeif (i + yBlockSize) < rows:ySize = yBlockSizeelse:ySize = rows - i# cycle col blocksizefor j in range(0, cols, xBlockSize):# set xblocksizeif (j + xBlockSize) < cols:xSize = xBlockSizeelse:xSize = cols - j# read block datadata = band.ReadAsArray(j, i, xSize, ySize).astype(np.float)sumIncludeZero = sumIncludeZero + np.sum(np.sum(data))mask = np.greater(data, 0)zeroCount = zeroCount + np.sum(np.sum(mask))# calculate
avg_includeZero = sumIncludeZero / (rows * cols)
avg_ignoreZero = sumIncludeZero / zeroCount
print('Average include zero: ', avg_includeZero)
print('Average ignore zero: ', avg_ignoreZero)

2）参考答案给的代码

# Assignment 4b:Write a script to calculate the average pixel value for the first band in aster.img
#reference answer
import os, sys
# import utilsuse_numeric = Truetry:from osgeo import ogr, gdalfrom osgeo.gdalconst import *import numpyos.chdir('E:/data/GDAL/ospy_data4')use_numeric = False
except ImportError:import ogr, gdalfrom gdalconst import *import Numericos.chdir('E:/data/GDAL/ospy_data4')# register all of the GDAL drivers
gdal.AllRegister()# open the image
ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:print('Could not open aster.img')sys.exit(1)# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount# # get the band and block sizes
# band = ds.GetRasterBand(1)
# blockSizes = utils.GetBlockSize(band)
# xBlockSize = blockSizes[0]
# yBlockSize = blockSizes[1]
band = ds.GetRasterBand(1)
xBlockSize = 20
yBlockSize = 20# initialize variables
count = 0
total = 0# loop through the rows
for i in range(0, rows, yBlockSize):if i + yBlockSize < rows:numRows = yBlockSizeelse:numRows = rows - i# loop through the columnsfor j in range(0, cols, xBlockSize):if j + xBlockSize < cols:numCols = xBlockSizeelse:numCols = cols - j# read the data and do the calculationsif use_numeric:data = band.ReadAsArray(j, i, numCols, numRows).astype(Numeric.Float)mask = Numeric.greater(data, 0)count = count + Numeric.sum(Numeric.sum(mask))total = total + Numeric.sum(Numeric.sum(data))else:data = band.ReadAsArray(j, i, numCols, numRows).astype(numpy.float)mask = numpy.greater(data, 0)count = count + numpy.sum(numpy.sum(mask))total = total + numpy.sum(numpy.sum(data))# print results
print('Ignoring 0:', total / count)
print('Including 0:', total / (rows * cols))

数据下载：https://download.csdn.net/download/qq_37935516/10795361

【GDAL学习】用GDAL读取栅格数据相关推荐

GDAL / OGR 学习手册 [02] ：栅格数据读取
目录一.栅格数据驱动二.gdal.Open 三.gdal.Dataset 四.获取影像的基本信息 1. 获取影像元数据 2. 获取影像基本信息一.栅格数据驱动 GDAL 通过数据驱动来识别各种类 ...
C# GDAL 学习一
最近一直琢磨如何用C#+GDAL读取栅格数据(.tif或.img),运气不错的在GDAL 的官网上找到一部分源码.经过本人测试,效果还不错.学习还将继续深入下去. 参考网址:http://trac.o ...
GDAL学习笔记02：GDAL基础知识
你的习惯决定了你会成为什么样的人. GDAL学习笔记02:GDAL基础知识前言 1. 版本 2. 摘要 3. 说明 4. 微信公众号GISRSGeography 一.GDAL简介二.导入GDAL ...
Python GDAL学习笔记(一)
目录一.读取tif格式影像导入GDAL模块并查看版本/位置打开影像查看行列号和波段数查看影像文件描述和其他属性打开数据集波段并获取统计信息利用numpy将波段数据写入数组将波段数据写入 ...
gdal 压缩tif_Python | GDAL处理影像
GDAL栅格数据处理栅格数据介绍栅格数据读取读取部分数据集坐标变换重采样什么是栅格数据基本上是一个大的二维或三维数组没有独立的几何对象,只有像素的集合二维:黑白图片三维:彩色/假彩 ...
R语言入门——读取栅格数据参数解读
读取栅格数据引言 1.数据读入 2.参数解读 2.1 class (类) 2.2 dimensions(维度) 2.3 resolution (像素) 2.4 extent(范围) 2.5 name ...
TF学习——TF数据读取：TensorFlow中数据读这三张图片的5个epoch +把读取的结果重新存到read 文件夹中
TF学习--TF数据读取:TensorFlow中数据读这三张图片的5个epoch +把读取的结果重新存到read 文件夹中目录实验展示代码实现实验展示代码实现 1.如果设置shuffle为T ...
ffplay.c学习-2-数据读取线程
ffplay.c学习-2-数据读取线程目录准备⼯作 avformat_alloc_context 创建上下⽂ ic->interrupt_callback avformat_open_inp ...
Adam学习25之读取sam生成的alignmentRecord含recordGroupDictionary
Adam学习25之读取sam生成的alignmentRecord含recordGroupDictionary 1.代码: package org.bdgenomics.adam.testimport ...

【GDAL学习】用GDAL读取栅格数据

【GDAL学习】用GDAL读取栅格数据相关推荐

最新文章

热门文章