【GIS】GDAL Python 影像裁剪
# -*- coding: utf-8 -*- """ Created on Fri Nov 30 11:45:03 2018 @author: Administrator """ from osgeo import gdal from osgeo import osr import numpy as np import math import time lonMeter = 0.00001141 latMeter = 0.00000899 #MeterParam = 0.00001 * 42496 / (124.44282531738276-124.3288421630859) MeterParam = 3.7282702222226876 def getSRSPair(dataset): ''' 获得给定数据的投影参考系和地理参考系 :param dataset: GDAL地理数据 :return: 投影参考系和地理参考系 ''' prosrs = osr.SpatialReference() prosrs.ImportFromWkt(dataset.GetProjection()) geosrs = prosrs.CloneGeogCS() return prosrs, geosrs def geo2lonlat(dataset, x, y): ''' 将投影坐标转为经纬度坐标(具体的投影坐标系由给定数据确定) :param dataset: GDAL地理数据 :param x: 投影坐标x :param y: 投影坐标y :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat) ''' prosrs, geosrs = getSRSPair(dataset) ct = osr.CoordinateTransformation(prosrs, geosrs) coords = ct.TransformPoint(x, y) return coords[:2] def lonlat2geo(dataset, lon, lat): ''' 将经纬度坐标转为投影坐标(具体的投影坐标系由给定数据确定) :param dataset: GDAL地理数据 :param lon: 地理坐标lon经度 :param lat: 地理坐标lat纬度 :return: 经纬度坐标(lon, lat)对应的投影坐标 ''' prosrs, geosrs = getSRSPair(dataset) ct = osr.CoordinateTransformation(geosrs, prosrs) coords = ct.TransformPoint(lon, lat) return coords[:2] def imagexy2geo(dataset, row, col): ''' 根据GDAL的六参数模型将影像图上坐标(行列号)转为投影坐标或地理坐标(根据具体数据的坐标系统转换) :param dataset: GDAL地理数据 :param row: 像素的行号 :param col: 像素的列号 :return: 行列号(row, col)对应的投影坐标或地理坐标(x, y) ''' trans = dataset.GetGeoTransform() px = trans[0] + col * trans[1] + row * trans[2] py = trans[3] + col * trans[4] + row * trans[5] return px, py def geo2imagexy(dataset, x, y): ''' 根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标(行列号) :param dataset: GDAL地理数据 :param x: 投影或地理坐标x :param y: 投影或地理坐标y :return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col) ''' trans = dataset.GetGeoTransform() a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]]) b = np.array([x - trans[0], y - trans[3]]) return np.linalg.solve(a, b) # 使用numpy的linalg.solve进行二元一次方程的求解 def imagexy2lonlat(dataset,row, col): ''' 影像行列转经纬度: :通过影像行列转平面坐标 :平面坐标转经纬度 ''' coords = imagexy2geo(dataset, row, col) coords2 = geo2lonlat(dataset,coords[0], coords[1]) return (coords2[0], coords2[1]) def lonlat2imagexy(dataset,x, y): ''' 影像行列转经纬度: :通过经纬度转平面坐标 :平面坐标转影像行列 ''' coords = lonlat2geo(dataset, x, y) coords2 = geo2imagexy(dataset,coords[0], coords[1]) return (int(round(abs(coords2[0]))), int(round(abs(coords2[1])))) if __name__ == '__main__': gdal.AllRegister() dataset = gdal.Open(r"D:\RSData\DAQING_SHAERTU\萨尔图区_大图:拼接\L19.tif") print('数据投影:') projection = dataset.GetProjection() print(projection) print('数据的大小(行,列):') print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize)) geotransform = dataset.GetGeoTransform() print('地理坐标:') print(geotransform) x = 464201 y = 5818760 lon = 122.47242 lat = 52.51778 row = 0 col = 0 # print('投影坐标 -> 经纬度:') # coords = geo2lonlat(dataset, x, y) # print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1])) # # print('经纬度 -> 投影坐标:') # coords = lonlat2geo(dataset, lon, lat) # print('(%s, %s)->(%s, %s)' % (lon, lat, coords[0], coords[1])) # # print('图上坐标 -> 投影坐标:') # coords = imagexy2geo(dataset, row, col) # print('(%s, %s)->(%s, %s)' % (row, col, coords[0], coords[1])) # # print('投影坐标 -> 图上坐标:') # coords = geo2imagexy(dataset, x, y) # print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1])) # print('图上坐标 -> 投影坐标:') # coords = imagexy2geo(dataset, row, col) # print('(%s, %s)->(%s, %s)' % (row, col, coords[0], coords[1])) # print('投影坐标 -> 经纬度:') # coords2 = geo2lonlat(dataset,coords[0], coords[1]) # print('(%s, %s)->(%s, %s)' % (coords[0], coords[1], coords2[0], coords2[1])) # coords = imagexy2lonlat(dataset, row, col) # print('影像像素 -> 经纬度:') # print('(%s, %s)->(%s, %s)' % ( row, col, coords[0], coords[1])) # coords = imagexy2lonlat(dataset, dataset.RasterXSize, dataset.RasterYSize) # print('影像像素 -> 经纬度:') # print('(%s, %s)->(%s, %s)' % ( dataset.RasterXSize, dataset.RasterYSize, coords[0], coords[1])) # # coords = lonlat2imagexy(dataset, 124.3288421630859, 46.391464001559044) # print('经纬度 -> 影像像素 :') # print('(%s, %s)->(%s, %s)' % ( 124.3288421630859, 46.391464001559044, coords[0], coords[1])) # coords = lonlat2imagexy(dataset, 124.44282531738276, 46.32796494040744) # print('经纬度 -> 影像像素 :') # print('(%s, %s)->(%s, %s)' % ( 124.44282531738276, 46.32796494040744, coords[0], coords[1])) #经纬度转像素 xoffset=0 yoffset=0 x,y = 125.059,46.894 xoffset,yoffset = lonlat2imagexy(dataset, x,y) print('坐标转换-对应行列像素位置') print('(%s, %s)->(%s, %s)' % (x,y, xoffset,yoffset)) width=int(500 * MeterParam) height=int(500 * MeterParam) if xoffset - width <= 0 and yoffset - height <= 0 : print("左上角") xoffset = 0 yoffset = 0 elif xoffset - width <= 0 and yoffset - height > 0 : print("左边") xoffset = 0 elif xoffset - width > 0 and yoffset - height <= 0 : print("顶边") yoffset = 0 else : print("中间区域") xoffset = xoffset - width yoffset = yoffset - height width = width * 2 height = height * 2 print('切割范围') print('宽高(%s, %s)->偏移起点(%s, %s)' % (width, height, xoffset,yoffset)) # xoffset,yoffset,width,height = 175360/2,123136/2,1000,1000 newData = np.zeros([width,height,3]) band = dataset.GetRasterBand(1) r=band.ReadAsArray(xoffset,yoffset,width,height) NoData = band.GetNoDataValue() newData[:,:,0] = r band = dataset.GetRasterBand(2) g=band.ReadAsArray(xoffset,yoffset,width,height) band = dataset.GetRasterBand(3) b=band.ReadAsArray(xoffset,yoffset,width,height) ticks = time.time() resultPath = "D:\\RS%s.jpg" % (int(ticks)) newData[:,:,0] = r newData[:,:,1] = g newData[:,:,2] = b format = "GTiff" driver = gdal.GetDriverByName(format) ds = driver.Create(resultPath, width, height, 3, gdal.GDT_Float32) geotransform1 = geotransform px = geotransform[0] + xoffset * geotransform[1] + yoffset * geotransform[2] py = geotransform[3] + xoffset * geotransform[4] + yoffset * geotransform[5] geotransform1 = (px, 0.29858214173896974, 0.0, py, 0.0, -0.29858214173896974) # print(geotransform1[0]) ds.SetGeoTransform(geotransform1) ds.SetProjection(projection) lay01= ds.GetRasterBand(1) lay02= ds.GetRasterBand(2) lay03= ds.GetRasterBand(3) # ds.GetRasterBand(1).SetNoDataValue(0) # ds.GetRasterBand(2).SetNoDataValue(0) # ds.GetRasterBand(3).SetNoDataValue(0) lay01.WriteArray(b) lay02.WriteArray(g) lay03.WriteArray(r) # ds.FlushCache() # ds = None del ds import cv2
import matplotlib.pyplot as plt img2=cv2.merge([r,g,b]) plt.imshow(img2) plt.xticks([]),plt.yticks([]) # 不显示坐标轴 plt.show() ticks = time.time() # cv2.imwrite("D:\\RS%s.jpg" % (int(ticks)) , img2) print("OK")
博客地址: http://www.cnblogs.com/defineconst/
博客版权: 本文以学习、研究和分享为主,欢迎转载和各类爬虫,但必须在文章页面明显位置给出原文链接。 如果文中有不妥或者错误的地方还望高手的您指出,以免误人子弟。如果您有更好的建议,不如留言一起讨论,共同进步! 再次感谢您耐心的读完本篇文章。
博客版权: 本文以学习、研究和分享为主,欢迎转载和各类爬虫,但必须在文章页面明显位置给出原文链接。 如果文中有不妥或者错误的地方还望高手的您指出,以免误人子弟。如果您有更好的建议,不如留言一起讨论,共同进步! 再次感谢您耐心的读完本篇文章。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· Linux系列:如何用 C#调用 C方法造成内存泄露
· AI与.NET技术实操系列(二):开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
· 探究高空视频全景AR技术的实现原理
· 理解Rust引用及其生命周期标识(上)
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· .NET 10首个预览版发布:重大改进与新特性概览!
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库