Python将普通图像转化为栅格影像

最新推荐文章于 2024-05-09 20:43:47 发布
最新推荐文章于 2024-05-09 20:43:47 发布
阅读量1.1k 收藏 19
点赞数 21
分类专栏： 遥感相关 python 文章标签： python arcgis
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_41099712/article/details/135300776
版权

引言

本人研究的方向是遥感，研究了2年也搞没清楚普通图像和遥感影像的区别，只知道到了多了地理坐标信息，但是经纬度信息映射到每个图像像素点的底层逻辑我还不太理解。因为现在需要使用python将图像转化为栅格影像，所以在此仔细研究一下。

地理坐标系和投影坐标系

【精选】地理坐标系VS投影坐标系_cgcs2000_3_degree_gk_cm_105e和gcs_china_geodetic_co ...

做遥感首先要理解地理坐标系（GeographicCoordinate System）和投影坐标系（ProjectedCoordinate System）2个概念。

  • 地理坐标系是球面坐标，参考平面是椭球面，坐标单位是经纬度；

  • 投影坐标系是平面坐标系，参考平面是水平面，坐标单位是米、千米等。

  • 地理坐标系转换到投影坐标系的过程理解为投影，即将不规则的地球曲面转换为平面。

img

以一个具体示例来初识ArcGIS中的坐标系，其全部参数拷贝在下面。这一示例是一个“投影坐标系（Projected Coordinate System）”。其名称是“WGS_1984_UTM_Zone_50N”。“WKID”是坐标系的编号，“ESPG”是“European Petroleum Survey Group”的缩写，表示其由“欧洲石油调查组织”发布。

可知，“WGS_1984_UTM_Zone_50N”这个投影坐标系由两部分组成：名为“Transverse_Mercator”的“投影（Projection）”和名为“GCS_WGS_1984”的“地理坐标系（GeographicCoordinate System）”。

WGS_1984_UTM_Zone_50N
WKID: 32650Authority: EPSG
Projection:Transverse_Mercator
False_Easting:500000.0
False_Northing:0.0
Central_Meridian:117.0
Scale_Factor:0.9996
Latitude_Of_Origin:0.0
Linear Unit: Meter(1.0)
GeographicCoordinate System: GCS_WGS_1984
Angular Unit:Degree (0.0174532925199433)
Prime Meridian:Greenwich (0.0)
Datum: D_WGS_1984
Spheroid: WGS_1984
Semimajor Axis: 6378137.0
Semiminor Axis: 6356752.314245179
Inverse Flattening: 298.257223563

地理坐标系

地理坐标系由三个参数来定义：角度单位（Angular Unit）、本初子午线（Prime Meridian）和大地测量系统（Datum）。

地理坐标系“GCS_WGS_1984”使用的角度单位为“度（Degree）”，0.0174532925199433这个数字等于“π/180”，使用的本初子午线为0.0度经线，即格林威治皇家天文台（Greenwich）所在位置的经线，使用的大地测量系统则为“D_WGS_1984”。
**地理坐标系的最重要的参数是“大地测量系统（Datum）”，而大地测量系统的最重要的参数是“椭球（Spheroid）”。**椭球相同，大地测量系统不一定相同，因为原点（origin）和方位（orientation）可以不同。想象一下，同一个椭球，首先可以固定在三维空间中的任意一个点，并且在固定于某点后还能以三个自由度任意地旋转其方位（朝向）。

WGS_1984”椭球的“长半轴（Semimajor Axis）”和“短半轴（Semiminor Axis）”分别为6378137.0和6356752.314245179，其“反扁率（Inverse Flattening）”为298.257223563，等于Semimajor Axis/( Semimajor Axis - Semiminor Axis)。

投影坐标系“WGS_1984_UTM_Zone_50N”使用的“投影（Projection）”名为“横轴墨卡托（Transverse_Mercator）”，然而这个名称并不能完全准确概括其投影。

投影坐标系

事实上，投影坐标系“WGS_1984_UTM_Zone_50N”这个名称中的“WGS_1984”指出了其地理坐标系为“GCS_WGS_1984”，而“UTM_Zone_50N”则指出了其投影。

“UTM_Zone_50N”这个名称指出，其投影方法是“通用横轴墨卡托（Universal Transverse Mercator，UTM）”，其投影带为北半球第50带，这个“Zone_50N”的“中央经线（Central Meridian）”正是117.0度，在“Transverse_Mercator”的参数中得到了体现。举一反三，“Xian_1980_GK_CM_117E”这个坐标系使用的地理坐标系为“GCS_Xian_1980”，而投影名称“GK_CM_117E”指出其使用以东经117度为中央经线的“高斯-克吕格（Gauss-Kruger，GK）”投影。投影的另一个重要参数是“东偏（False Easting）”。

有些投影会在X坐标值前加上投影带号，比如：“Xian_1980_GK_Zone_20”的“false_easting”参数为20500000.0，其中20为投影带号，而“Xian_1980_GK_CM_117E”的“false_easting”参数为500000.0，尽管它们的中央经线都为东经117度。

地理坐标系和投影坐标区别及常用操作: https://zhuanlan.zhihu.com/p/484585670

GeoTIFF数据格式

GeoTIFF是一种常用的地理信息系统（GIS）数据格式，它结合了TIFF（Tagged Image File Format）图像格式和地理空间信息。GeoTIFF文件可以包含地理坐标系、投影信息、地理位置和其他与地理空间数据相关的元数据。

GeoTIFF文件可以存储各种类型的地理空间数据，包括栅格数据（如卫星影像、数字高程模型等）和矢量数据（如点、线、面等地理要素）。它支持不同的地理坐标系和投影系统，可以准确地表示地球上的位置和形状。

GeoTIFF文件的优点包括：

  1. 与TIFF格式兼容：GeoTIFF文件可以使用标准的TIFF查看器和编辑器打开和编辑。
  2. 地理空间信息：GeoTIFF文件可以包含地理坐标系、投影信息和其他与地理空间数据相关的元数据，使得数据可以准确地在地球上的位置上显示和分析。
  3. 多波段支持：GeoTIFF文件可以存储多个波段的数据，例如多光谱卫星影像，方便进行多光谱分析。
  4. 跨平台兼容性：GeoTIFF文件可以在不同的GIS软件和平台上使用，包括ArcGIS、QGIS、ENVI等。

使用GeoTIFF文件，可以进行地理空间数据的获取、处理、分析和可视化，是GIS领域中常用的数据格式之一。

以下是一些常见的地理空间信息在GeoTIFF文件中的保存方式：

  1. 坐标系信息：GeoTIFF文件中保存了坐标系的相关信息，包括坐标系的名称、投影信息、椭球体参数等。这些信息可以通过GeoTIFF文件的元数据标签来获取。
  2. 地理范围：GeoTIFF文件可以保存地理范围的信息，即图像所覆盖的地理区域的边界。通常使用左上角和右下角的经纬度坐标来表示地理范围。
  3. 分辨率：GeoTIFF文件可以保存图像的分辨率信息，即每个像素代表的地理距离。通常以米或度为单位。
  4. 像素对应的地理坐标：GeoTIFF文件可以保存每个像素对应的地理坐标信息。这样可以通过像素坐标来获取对应的地理位置。

这些地理空间信息是通过GeoTIFF文件的元数据标签来保存的，常见的元数据标签包括"GeoKeyDirectoryTag"、“ModelPixelScaleTag”、"ModelTiepointTag"等。可以使用专门的GIS软件或者Python库（如GDAL）来读取和处理GeoTIFF文件的地理空间信息。

GDAL库

GDAL 官方文档 (osgeo.cn)

GDAL（Geospatial Data Abstraction Library）是一个开源的地理空间数据处理库，它提供了对各种格式的地理空间数据进行读取、写入和转换的功能。GDAL库支持包括GeoTIFF、Shapefile、NetCDF、HDF等在内的多种地理空间数据格式。

GDAL库的主要功能包括：

  1. 数据读取和写入：GDAL可以读取和写入各种地理空间数据格式，包括栅格数据和矢量数据。它可以从文件、数据库或内存中读取数据，并将数据写入到文件或其他数据源中。

  2. 数据转换和投影：GDAL可以进行地理空间数据的投影转换、坐标系转换和数据格式转换。它提供了丰富的投影和坐标系操作功能，可以将数据从一个坐标系转换到另一个坐标系，并进行投影变换。

  3. 数据操作和分析：GDAL提供了一些基本的数据操作和分析功能，如裁剪、重采样、统计等。它可以对栅格数据进行像素级别的操作，如计算均值、最大值、最小值等。

  4. 数据集信息获取：GDAL可以获取地理空间数据集的元数据信息，包括坐标系信息、地理范围、分辨率等。它还可以获取数据集中的波段信息、属性表信息等。

GDAL是一个功能强大且广泛应用的地理空间数据处理库，它提供了丰富的功能和灵活的接口，适用于各种地理空间数据处理需求。在Python中，可以使用GDAL库的Python绑定（即GDAL Python）来进行地理空间数据处理。

基础操作

from osgeo import gdal

# 打开影像文件
dataset = gdal.Open('data/other/guoyang2019.tif')
# 获取影像宽度和高度
width = dataset.RasterXSize
height = dataset.RasterYSize
# 获取影像波段数
band_count = dataset.RasterCount
# 获取投影信息
projection = dataset.GetProjection()
# 获取地理转换信息
geotransform = dataset.GetGeoTransform()
# 打印影像基本信息
print("影像宽度:", width)
print("影像高度:", height)
print("波段数:", band_count)
print("投影信息:", projection)
print("地理转换信息:", geotransform)
# 关闭影像文件
dataset = None

输出：

影像宽度: 7471
影像高度: 4896
波段数: 1
投影信息: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
地理转换信息: (115.878317255137, 8.982254799999953e-05, 0.0, 33.785346925021, 0.0, -8.982254799999988e-05)

其中地理转换信息即仿射矩阵，实现了从像素点到经纬度坐标的变换，其含义为

GeoTransform = [115.878317255137, 8.982254799999953e-05, 0.0, 33.785346925021, 0.0, -8.982254799999988e-05]
'''
0：图像左上角的X坐标；
1：图像东西方向分辨率；
2：旋转角度，如果图像北方朝上，该值为0；
3：图像左上角的Y坐标；
4：旋转角度，如果图像北方朝上，该值为0；
5：图像南北方向分辨率；
'''
得到这六个参数之后就可以进行图像像素坐标（即行列号）和地理坐标之间的变换：
#像素坐标和地理坐标仿射变换
XGeo = GeoTransform[0]+GeoTransform[1]*Xpixel+Ypixel*GeoTransform[2];
YGeo = GeoTransform[3]+GeoTransform[4]*Xpixel+Ypixel*GeoTransform[5];

读取和写入文件的函数：

def read_img(filename):
        dataset = gdal.Open(filename)
        im_width = dataset.RasterXSize #栅格矩阵的列数
        im_height = dataset.RasterYSize #栅格矩阵的行数
        im_geotrans = dataset.GetGeoTransform() #仿射矩阵
        im_proj = dataset.GetProjection() #地图投影
        im_data = dataset.ReadAsArray(0,0,im_width,im_height)
        del dataset
        return im_proj,im_geotrans,im_data
def write_img(filename,im_proj,im_geotrans,im_data):
        if "int8" in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        if "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

        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

实战

写这篇文章的目的是因为我目前深度学习拼接上的作物分布图只是一张PNG图片，没有地理坐标信息，现在需要通过代码添加上：

SoyaMap

from osgeo import gdal, osr

path = f"img/predict2020/SoyaMap.png"
# 读取图像
image = Image.open(path)
# 转换为numpy数组
image_array = (1-np.array(image)/255).astype(np.int8)
# 创建一个空间参考对象 ,将空间坐标系设置为EPSG:4326
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
# 地理转换信息,图像左上角的经纬度坐标点和像素大小
geotransform = [115.86991557438377,8.983152841195215E-5, 0.0, 33.789500591456914, 0.0, -8.983152841195215E-5]
write_img("img/predict2020/SoyaMap.tif",srs.ExportToWkt(),geotransform,image_array)

导出成功后，使用Arcgis打开制图，简直不要太酷啦。

涡阳大豆分布图2020

齐硕君
关注
  • 21
    点赞
  • 踩
  • 19
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
python+GDAL实现将GeoJson转换为栅格影像
做GIS梦的人
08-21 481
最近遇到一个需求，需要将GeoJson转为jpg图片或者tiff格式的二值图像。这个需求，在web前端实现肯定是没有问题的，可以基于WebGIS框架进行渲染GeoJson，然后canvas截图就可以。不过对于一个后端项目来说，这样的实现方式就稍微有些麻烦了，效率也低。调研了一些网上的方法，大部分是用来实现GeoJson转栅格的，但是看到他们的代码里，总要先设置一个模版栅格文件，然后从模版文件获取像素数量、坐标系、transform参数。这样的方法对于生成固定分辨率和像元数量的需求是可以满足的。
关于Python 的简单栅格图像边界提取方法
09-19
今天小编就为大家分享一篇关于Python 的简单栅格图像边界提取方法，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧
利用Python栅格化地图，含代码（以成都市为例）
weixin_46357639的博客
11-01 1541
利用python将研究区域进行栅格化
【Python&RS】基于Python栅格数据/遥感影像投影转换
Remote Sensing迷途小书童
06-28 799
为了解放双手方便批量投影转换，最近同步研究了一下如何利用Python实现遥感影像的投影转换，主打的就是一个懒。博文基于GDAL实现了栅格数据的投影转换！
Python批量实现栅格影像可视化
qq_43874102的博客
04-22 3293
任务描述： 对栅格影像进行可视化并保存图片 数据准备： 文件夹(工作空间)内含有需要批量可视化的tif影像 工具准备： Python环境 操作： 替换代码中的路径、根据需要调整可视化图片格式 # -*- coding: utf-8 -*- import os import numpy as np from osgeo import gdal import earthpy_plot_revised as ep # 文件目录 root_path = r"D:\DATA\04_T" out_path = r"
Python Gdal 栅格数据处理之将数组存为Tif数据
weixin_43896283的博客
06-21 2016
Python Gdal 栅格数据处理之将数组存为Tif数据
python gdal 选择通道将tif格式影像转换为png
qq_43509463的博客
11-02 2905
def RGBNir_tiftoRGB_PNG(): options = gdal.TranslateOptions(format='PNG', bandList=[4, 3, 2]) gdal.Translate('RGB.png', 'GF.tif', options=options) 其中bandlist为原tif中用于生成png的波段数组，如缺省该参数，则默认为[1, 2, 3] 以GF-2影像为例，选择波段[4 , 3, 2]，生成的PNG如下 如由其他要求可自行前往源代.
Python遥感图像处理应用篇(十六)：GDAL 将归一化处理csv数据转化为遥感影像
soderayer的博客
06-16 1262
将上一篇文章中得到的计算结果作为转换数据。链接如下：Python遥感图像处理应用篇(十五)：GDAL 读取多光谱数据做归一化处理_空中旋转篮球的博客-CSDN博客基本思路，先读取csv数据，现读取每一列数据为一维数组，然后根据原始影像数据的大小（row，col）将一维数组塑形，即使用reshape函数变换为（row，col）形状对的二维数组。使用gdal将二维数组转换为tif图像。完整代码如下： 修改一下以上代码，可以将数据转化为和原始数据相同波段数的多波段遥感影像。可以使用for循环一个通道一个通道的
【Python&RS】GDAL批量裁剪遥感影像/栅格数据
Remote Sensing迷途小书童
05-26 1959
GDAL批量裁剪遥感影像/栅格数据
Python中ArcPy实现栅格图像文件由HDF格式批量转换为TIFF格式
疯狂学习GIS的博客
04-15 4504
本文介绍基于Python中Arcpy模块，实现HDF格式栅格图像文件批量转换为TIFF格式的方法~
Python叠加两幅栅格图像的实现方法
09-19
今天小编就为大家分享一篇Python叠加两幅栅格图像的实现方法，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧
python2.7栅格数据批量转换投影
04-29
python2.7栅格数据批量转换投影:ProjectRaster_management (in_raster, out_raster, out_coor_system, {resampling_type}, {cell_size}, {geographic_transform}, {Registration_Point}, {in_coor_system}) ...
vectorize:将栅格图像转换为矢量化轮廓或多边形
05-13
向量化将栅格图像转换为矢量化轮廓或多边形
【 书生·浦语大模型实战营】作业（六）：Lagent & AgentLego 智能体应用搭建
专注大数据与人工智能技术分享，欢迎私信加群互相学习！
05-03 1549
本篇笔记内容主要为 【书生·浦语大模型实战营】作业（六）：Lagent & AgentLego 智能体应用搭建智能体的作业 ，主要对智能体应用理论部分及具体搭建流程详细介绍，并分别进行实战应用及可视化展示，欢迎大家交流学习！
第五节课《LMDeploy 量化部署 LLM 实践》
fanre的专栏
05-06 630
PDF链接：https://pan.baidu.com/s/1JFtvBWgEGFWJq8pHafvIUg?pwd=6666提取码：6666。
Python从0到POC编写--实用小脚本
最新发布
余十步的博客
05-09 920
模块中定义了一些常用的属性（包含所有数字，字母，可打印的所有ascii码等）第二步就是要把这些没有写的协议，或者没有补全的协议全部替换为我们的内容，在ctf比赛中，为了能够快速获取flag，获得更多的分数，现在呢，代码是这样子的，需要把一些标签替换为空，把。这份列表呢，奇奇怪怪，有些写错了，有些没补全。不管用户输入的是什么，都替换为我们安排的内容，现在呢，我们有了一份奇怪的列表要去判断，我们只知道前三位，然后需要去爆破后两位，当然，比赛一般有两种获取flag的方式，在解题的过程中，只能看到一半的密码，
继承与组合
g1446992400的博客
05-06 949
继承就是创建新类的一种方式，这个新类可以继承一个或者多个其他类的属性新的类如果有自己属性，那就叫派生如果属性不封的情况下，谁实例化得到的self 就去谁里面找如果属性封装的情况下 ，谁实例化得到的self 无效，只能在当前所在的类的名称空间里面找。
十二届蓝桥杯Python组1月中/高级试题 第四题
走出那座山
05-09 468
输入一个正整数n，如果n为偶数，则计算1/2+1/4+1/6…+1/n的结果，如果n为奇数时，则计 算1/1+1/3+1/5…输出最后计算的结果。输入描述：输入一个正整数n 输出描述：输出最后计算的结果（最后结果保留两位小数，不要在每次除法时保留小数）输入为5，则计算1/1+1/3+1/5的结果，输出为1.53。例如：输入为4，则计算1/2+1/4的结果，输出为0.75。第四题（难度系数 4，30 个计分点）
python栅格图像
03-20
栅格图像是由像素组成的二维数组，每个像素都有一个特定的数值表示其属性或特征。Python提供了多种库和工具来处理栅格图像，例如NumPy、PIL（Python Imaging Library）和OpenCV等。下面是一个使用NumPy库处理栅格图像的例子[^2]： ```python import numpy as np # 创建一个3x3的栅格图像 image = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) # 获取图像的形状（行数和列数） rows, cols = image.shape # 访问图像的像素值 pixel_value = image[1, 2] # 获取第2行第3列的像素值（从0开始计数） # 修改图像的像素值 image[0, 0] = 10 # 将第1行第1列的像素值修改为10 # 遍历图像的所有像素 for i in range(rows): for j in range(cols): print(image[i, j]) # 对图像进行统计分析 mean = np.mean(image) # 计算图像的平均值 std = np.std(image) # 计算图像的标准差 max_value = np.max(image) # 计算图像的最大值 min_value = np.min(image) # 计算图像的最小值 ``` 这是一个简单的栅格图像处理的例子，你可以根据具体的需求使用不同的库和方法来处理栅格图像。如果你有其他关于栅格图像的问题，请告诉我。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值