python+GDAL+numpy，点图层提取栅格像元数据

这部强调：投影坐标一定要一致（shp和栅格）！！！投影坐标一定要一致（shp和栅格）！！！投影坐标一定要一致（shp和栅格）！！！CRS.from_epsg('32650')!CRS.from_epsg('32650')!!CRS.from_epsg('32650')!!

• EPSG:32650: WGS 84 / UTM zone 50N

好了继续，有几个办法，一个是用gdal readRaster，或者把栅格转数组。。。读对应位置的数据（注意位置要对应上）

from osgeo import gdal,ogr

import struct

src_filename = 'D:/Thesis/ML/aodrepro/MCD19A2.A2018001.h28v06.006.2018121012322.hdf.tif'

shp_filename = 'D:/Thesis/point/point72.shp'

src_ds=gdal.Open(src_filename)

gt=src_ds.GetGeoTransform()

rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)

lyr=ds.GetLayer()

plist=list()

#inyy=[]

for feat in lyr:

geom = feat.GetGeometryRef()

mx,my=geom.GetX(), geom.GetY() #coord in map units

#Convert from map to pixel coordinates.

#Only works for geotransforms with no rotation.

px = int((mx - gt[0]) / gt[1]) #x pixel

py = int((my - gt[3]) / gt[5]) #y pixel----- ##实在不行就用数组提取吧

# band = src_ds.GetRasterBand(1)#2

# arr = band.ReadAsArray()

# inyy=arr[py,px]

# plist.append(inyy)

# print (inyy)

structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 64 bit int aka 'double'px和py是从左到右，从下到上逐个计算位置的个数，其源码是int，表示个数

intval = struct.unpack('h' , structval) #use the 'double' format code 4 bytes

plist.append(intval[0])

###structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) 解释一下，px是算的，见上面公式，是坐标减去栅格最左值，除以像元大小，就是第几个像元了，同理，py；1，1是计算一个像元的意思，横着1，竖着1.。。。后面就是16位，，，，

得到的结果可以存到csv里

：

1 import numpy as np
2 np.savetxt('E:\\forpython\\featvector.csv',data_to_save,delimiter=',')