Geoprocessing with Python （一）

最新推荐文章于 2023-03-17 08:46:22 发布

wblong_cs

最新推荐文章于 2023-03-17 08:46:22 发布

阅读量790

点赞数

分类专栏： python_gis

本文链接：https://blog.csdn.net/mrbaolong/article/details/104241068

版权

python_gis 专栏收录该内容

20 篇文章 7 订阅

订阅专栏

Week 7: Misc stuff

Convert a shapefile to KML

Use “STATE_NAME” as the name field for the KML features
Include all of the attributes in the description

第七周：杂项

shapefile 转换 KML

使用“STATE_NAME”作为 KML要素的Name字段
在描述中包含所有的属性

ogr2ogr -f "KML" -dsco NameField="STATE_NAME" states.kml states.shp

Week 6 : More raster processing

Use a 3x3 high-pass filter to detect edges in band 1 of smallaster.img.

Use pixel notation (this is why you’re using smallaster.img – it would take forever to run on aster.img).
The output type will be float.
Don’t build pyramids so you can easily compare your output with output.img.

Use a 3x3 high-pass filter to detect edges in band 1 of aster.img.

Use slice notation. It will take a minute or two to run, so you should probably test on smallaster.img and then run in on aster.img after everything seems to be working.
The output type will be float.
Don’t build pyramids so you can easily compare your output with output.img.

更多栅格影像处理

使用 3x3 高通滤波来检测 smallaster.img （波段1）边缘。

采用像素表示法。
输出类型 float。
不建立金子塔影像

使用 3x3 高通滤波来检测 aster.img （波段1）边缘。

采用切片表示法。
输出类型 float。
不建立金子塔影像

# homework 6a
# high-pass filter using pixel notation
import os, sys, time
import numpy as np

from osgeo import gdal
from osgeo.gdalconst import *

import matplotlib.pyplot as plt
#plt.rcParams['figure.dpi'] = 300 #分辨率

start = time.time()

#os.chdir('c:/temp/py6')
# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open('smallaster.img', GA_ReadOnly)
if inDs is None:
  print ('Could not open smallaster.img')
  sys.exit(1)

# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize

# read the input data
inBand = inDs.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, cols, rows).astype(np.float32)

# do the calculation
outData = np.zeros((rows, cols), np.float32)

for i in range(1, rows-1):
    for j in range(1, cols-1):
        outData[i,j] = ((-0.7 * inData[i-1,j-1]) + 
                        (-1.0 * inData[i-1,j])   + 
                        (-0.7 * inData[i-1,j+1]) + 
                        (-1.0 * inData[i,j-1])   + 
                        ( 6.8 * inData[i,j])     + 
                        (-1.0 * inData[i,j+1])   +
                        (-0.7 * inData[i+1,j-1]) + 
                        (-1.0 * inData[i+1,j])   + 
                        (-0.7 * inData[i+1,j+1]))
# create the output image
driver = inDs.GetDriver()
outDs = driver.Create('highpass1.img', cols, rows, 1, GDT_Float32)
if outDs is None:
  print ('Could not create highpass1.img')
  sys.exit(1)
outBand = outDs.GetRasterBand(1)

# write the output data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
stats = outBand.GetStatistics(0, 1)

# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())

print (time.time() - start, 'seconds')

plt.subplot(121)
plt.imshow(inData,cmap='gray')
plt.subplot(122)
plt.imshow(outData,cmap='gray')
plt.savefig("res1.png")
inDs = None
outDs = None

74.08075737953186 seconds

在这里插入图片描述

# homework 6b
# high-pass filter using slice notation
import os, sys, time
import numpy as np

from osgeo import gdal
from osgeo.gdalconst import *

import matplotlib.pyplot as plt

#plt.rcParams['figure.dpi'] = 300 #分辨率

start = time.time()

#os.chdir('c:/temp/py6')

# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open('smallaster.img', GA_ReadOnly)
if inDs is None:
  print ('Could not open aster.img')
  sys.exit(1)

# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize

# read the input data
inBand = inDs.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, cols, rows).astype(np.uint32)


# do the calculation
outData = np.zeros((rows, cols), np.float32)

outData[1:rows-1,1:cols-1] = ((-0.7 * inData[0:rows-2,0:cols-2]) +
                              (-1.0 * inData[0:rows-2,1:cols-1]) +
                              (-0.7 * inData[0:rows-2,2:cols]) +
                              (-1.0 * inData[1:rows-1,0:cols-2]) + 
                              ( 6.8 * inData[1:rows-1,1:cols-1]) +
                              (-1.0 * inData[1:rows-1,2:cols]) + 
                              (-0.7 * inData[2:rows,0:cols-2]) +
                              (-1.0 * inData[2:rows,1:cols-1]) + 
                              (-0.7 * inData[2:rows,2:cols]))

# outData[1:rows-1,1:cols-1] = ((-1.0 * inData[0:rows-2,0:cols-2]) +
#                               (-1.0 * inData[0:rows-2,1:cols-1]) +
#                               (-1.0 * inData[0:rows-2,2:cols]) +
#                               (-1.0 * inData[1:rows-1,0:cols-2]) + 
#                               ( 8.0 * inData[1:rows-1,1:cols-1]) +
#                               (-1.0 * inData[1:rows-1,2:cols]) + 
#                               (-1.0 * inData[2:rows,0:cols-2]) +
#                               (-1.0 * inData[2:rows,1:cols-1]) + 
#                               (-1.0 * inData[2:rows,2:cols]))


# create the output image
driver = inDs.GetDriver()
outDs = driver.Create('highpass3.img', cols, rows, 1, GDT_Float32)#GDT_CInt32
if outDs is None:
  print ('Could not create highpass3.img')
  sys.exit(1)
outBand = outDs.GetRasterBand(1)

# write the output data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
stats = outBand.GetStatistics(0, 1)

# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())

print (time.time() - start, 'seconds')

plt.subplot(121)
plt.imshow(inData,cmap='gray')

plt.subplot(122)
plt.imshow(outData,cmap='gray')
plt.savefig("res2.png")
inDs = None
outDs = None

0.48896169662475586 seconds

在这里插入图片描述

# homework 6b
# low-pass filter using slice notation
import os, sys, time
import numpy as np

from osgeo import gdal
from osgeo.gdalconst import *

import matplotlib.pyplot as plt

#plt.rcParams['figure.dpi'] = 300 #分辨率

start = time.time()

#os.chdir('c:/temp/py6')

# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open('smallaster.img', GA_ReadOnly)
if inDs is None:
  print ('Could not open aster.img')
  sys.exit(1)

# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize

# read the input data
inBand = inDs.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, cols, rows).astype(np.uint32)

#plt.figure(figsize=(10,10)) 


# do the calculation
outData = np.zeros((rows, cols), np.float32)

outData[1:rows-1,1:cols-1] = ((0.111 * inData[0:rows-2,0:cols-2]) +
                              (0.111 * inData[0:rows-2,1:cols-1]) +
                              (0.111* inData[0:rows-2,2:cols]) +
                              (0.111 * inData[1:rows-1,0:cols-2]) + 
                              (0.111 * inData[1:rows-1,1:cols-1]) +
                              (0.111 * inData[1:rows-1,2:cols]) + 
                              (0.111 * inData[2:rows,0:cols-2]) +
                              (0.111 * inData[2:rows,1:cols-1]) + 
                              (0.111 * inData[2:rows,2:cols]))

# outData[1:rows-1,1:cols-1] = ((-1.0 * inData[0:rows-2,0:cols-2]) +
#                               (-1.0 * inData[0:rows-2,1:cols-1]) +
#                               (-1.0 * inData[0:rows-2,2:cols]) +
#                               (-1.0 * inData[1:rows-1,0:cols-2]) + 
#                               ( 8.0 * inData[1:rows-1,1:cols-1]) +
#                               (-1.0 * inData[1:rows-1,2:cols]) + 
#                               (-1.0 * inData[2:rows,0:cols-2]) +
#                               (-1.0 * inData[2:rows,1:cols-1]) + 
#                               (-1.0 * inData[2:rows,2:cols]))


# create the output image
driver = inDs.GetDriver()
outDs = driver.Create('lowpass3.img', cols, rows, 1, GDT_Float32)#GDT_CInt32
if outDs is None:
  print ('Could not create highpass3.img')
  sys.exit(1)
outBand = outDs.GetRasterBand(1)

# write the output data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
stats = outBand.GetStatistics(0, 1)

# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())



print (time.time() - start, 'seconds')

plt.subplot(121)
plt.imshow(inData,cmap='gray')

plt.subplot(122)
plt.imshow(outData,cmap='gray')

plt.savefig("res3.png")
inDs = None
outDs = None

0.4882187843322754 seconds

在这里插入图片描述

Projecting rasters

设置栅格投影

from osgeo import gdal, osr 
from osgeo.gdalconst import *
import time
import matplotlib.pyplot as plt 

start = time.time()

#获取源数据及栅格信息
gdal.AllRegister()
src_data = gdal.Open("smallaster.img")
#获取源的坐标信息
srcSRS_wkt=src_data.GetProjection()
srcSRS=osr.SpatialReference()
srcSRS.ImportFromWkt(srcSRS_wkt)
#print("原数据空间参考：",srcSRS)
#获取栅格尺寸
src_width = src_data.RasterXSize
src_height = src_data.RasterYSize
src_count=src_data.RasterCount
#print("原数据图像信息：\n",src_width,src_height,src_count)
#获取源图像的仿射变换参数
src_trans=src_data.GetGeoTransform()

#左上角空间坐标
OriginLX_src=src_trans[0]
OriginTY_src=src_trans[3]
# 像素比
pixl_w_src=src_trans[1]
pixl_h_src=src_trans[5]
# 右下角坐标
OriginRX_src=OriginLX_src+pixl_w_src*src_width
OriginBY_src=OriginTY_src+pixl_h_src*src_height
#创建输出图像
driver = gdal.GetDriverByName("GTiff")
driver.Register()
dst_data = driver.Create("tpix1.tif",src_width,src_height,src_count)

#设置输出图像的空间参考WGS84
dstSRS=osr.SpatialReference()
dstSRS.ImportFromEPSG(4326)
#投影转换
ct=osr.CoordinateTransformation(srcSRS,dstSRS)
#计算目标影像的左上和右下坐标,即目标影像的仿射变换参数
OriginLX_dst,OriginTY_dst,temp=ct.TransformPoint(OriginLX_src,OriginTY_src)
OriginRX_dst,OriginBY_dst,temp=ct.TransformPoint(OriginRX_src,OriginBY_src)
 
pixl_w_dst=(OriginRX_dst-OriginLX_dst)/src_width
pixl_h_dst=(OriginBY_dst-OriginTY_dst)/src_height

dst_trans=[OriginLX_dst,pixl_w_dst,0,OriginTY_dst,0,pixl_h_dst]
#
dstSRS_wkt=dstSRS.ExportToWkt()
#设置仿射变换系数及投影
dst_data.SetGeoTransform(dst_trans)
dst_data.SetProjection(dstSRS_wkt)
#重新投影
gdal.ReprojectImage(src_data,dst_data,srcSRS_wkt,dstSRS_wkt,GRA_Bilinear)
print (time.time() - start, 'seconds')
#创建四级金字塔
#gdal.SetConfigOption('HFA_USE_RRD', 'YES')
#dst_data.BuildOverviews(overviewlist=[2,4,8,16])
#计算统计值
# for i in range(0,src_count):
#     band=dst_data.GetRasterBand(i+1)
#     band.SetNoDataValue(-99)
#     print (band.GetStatistics(0,1))

# read the input data
inBand = src_data.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, src_width,src_height).astype(np.uint32)
plt.subplot(121)
plt.imshow(inData)

inBand = dst_data.GetRasterBand(1)
dstData = inBand.ReadAsArray(0, 0, src_width,src_height).astype(np.uint32)
plt.subplot(122)
plt.imshow(dstData)

1.7203800678253174 seconds

from osgeo import gdal, osr 
from osgeo.gdalconst import *
import numpy as np
import time

start = time.time()

#获取源数据及栅格信息
gdal.AllRegister()
src_data = gdal.Open("smallaster.img")
#获取源的坐标信息
srcSRS_wkt=src_data.GetProjection()
srcSRS=osr.SpatialReference()
srcSRS.ImportFromWkt(srcSRS_wkt)
#print("原数据空间参考：",srcSRS)
#获取栅格尺寸
src_width = src_data.RasterXSize
src_height = src_data.RasterYSize
src_count=src_data.RasterCount
#print("原数据图像信息：\n",src_width,src_height,src_count)
#获取源图像的仿射变换参数
src_trans=src_data.GetGeoTransform()

#左上角空间坐标
OriginLX_src=src_trans[0]
OriginTY_src=src_trans[3]
# 像素比
pixl_w_src=src_trans[1]
pixl_h_src=src_trans[5]
# 右下角坐标
OriginRX_src=OriginLX_src+pixl_w_src*src_width
OriginBY_src=OriginTY_src+pixl_h_src*src_height
#创建输出图像
driver = gdal.GetDriverByName("GTiff")
driver.Register()
dst_data = driver.Create("tpix2.tif",src_width,src_height,src_count)

for i in range(src_count):
    inBand = src_data.GetRasterBand(i+1)
    inData = inBand.ReadAsArray(0, 0, src_width, src_height).astype(np.uint32)
    outBand = dst_data.GetRasterBand(i+1)
    outBand.WriteArray(inData, 0, 0)# write the output data
    outBand.FlushCache()
    stats = outBand.GetStatistics(0, 1)
    inData=None

#设置输出图像的空间参考WGS84
dstSRS=osr.SpatialReference()
dstSRS.ImportFromEPSG(4326)
#投影转换
ct=osr.CoordinateTransformation(srcSRS,dstSRS)
#计算目标影像的左上和右下坐标,即目标影像的仿射变换参数
OriginLX_dst,OriginTY_dst,temp=ct.TransformPoint(OriginLX_src,OriginTY_src)
OriginRX_dst,OriginBY_dst,temp=ct.TransformPoint(OriginRX_src,OriginBY_src)
 
pixl_w_dst=(OriginRX_dst-OriginLX_dst)/src_width
pixl_h_dst=(OriginBY_dst-OriginTY_dst)/src_height

dst_trans=[OriginLX_dst,pixl_w_dst,0,OriginTY_dst,0,pixl_h_dst]
#
dstSRS_wkt=dstSRS.ExportToWkt()
#设置仿射变换系数及投影
dst_data.SetGeoTransform(dst_trans)
dst_data.SetProjection(dstSRS_wkt)
print (time.time() - start, 'seconds')
#重新投影
#gdal.ReprojectImage(src_data,dst_data,srcSRS_wkt,dstSRS_wkt,GRA_Bilinear)
#创建四级金字塔
#gdal.SetConfigOption('HFA_USE_RRD', 'YES')
#dst_data.BuildOverviews(overviewlist=[2,4,8,16])
#计算统计值
# for i in range(0,src_count):
#     band=dst_data.GetRasterBand(i+1)
#     band.SetNoDataValue(-99)
#     print (band.GetStatistics(0,1))

0.31161069869995117 seconds

Method comparison

比较运算

# homework 6b
# low-pass filter using slice notation
import os, time,sys
import numpy as np

from osgeo import gdal
from osgeo.gdalconst import *

import matplotlib.pyplot as plt

start = time.time()

# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open('dem.tiff', GA_ReadOnly)
if inDs is None:
  print ('Could not dem.tiff')
  sys.exit(1)

# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize
nbands = inDs.RasterCount
print(rows,cols,nbands)
# read the input data
inBand = inDs.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, cols, rows).astype(np.uint32)

plt.figure(figsize=(10,10))
plt.subplot(141)
plt.imshow(inData,cmap='gray')
plt.title("all")

outData = np.zeros((rows, cols), np.uint32)
outData1 = np.greater(inData, 1000)

plt.subplot(142)
plt.imshow(outData1,cmap='gray')
plt.title("elev>1000")

outData2= np.where((inData < 1000) & ( inData > 500), 1, 0)

plt.subplot(143)
plt.imshow(outData2,cmap='gray')
plt.title("elev>500 & elev <1000")

plt.savefig("res4.png")

324 368 1

在这里插入图片描述
参考文档：

Geoprocessing with Python using Open Source GIS(https://www.gis.usu.edu/~chrisg/python/2009/)

wblong_cs

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
0
评论
Geoprocessing with Python （一）

Week 7: Misc stuffConvert a shapefile to KMLUse “STATE_NAME” as the name field for the KML featuresInclude all of the attributes in the description第七周：杂项shapefile 转换 KML使用“STATE_NAME”作为 KML要素...
复制链接

扫一扫