SAR 图像像素点转化经纬度

有的 tiff 格式的 SAR 图像其可能存在 CRS 和转化矩阵，可以使用这两个参数将 SAR 图像中的像素点转化为实际的经纬度。

CRS 地理坐标是用来表示地球上的具体某一点的坐标系，类似于二维平面下有极坐标系和直角坐标系一样，CRS 也有很多种不同的坐标系，下面输出的 CRS 就是数据集使用的地理坐标系，但是这个不是标准的经纬度坐标，输出 CRS 信息是用来描述这个坐标系的具体参数，但是实际使用的时候并不需要全部了解这些参数的意义，直接整体使用即可。

将像素点转化为经纬度有两步仿射变化：

  1. 使用 transform 矩阵将像素点转化为该数据集下的 CRS 地理坐标系。
  2. 将该数据集的坐标系转化为标准经纬度坐标系。

第一步的 transform 矩阵是 tiff 图像种可以直接获取的，如下面所示。但是如果输出的是单位矩阵的话表示该数据集没有可使用的转化矩阵，没办法转化为 CRS 坐标。

1
2
3
4
5
6
7
8
9
10
11
import rasterio

# 打开 TIFF 文件
with rasterio.open('../dataset/AIR-SARShip-1.0/data_spilt/images/SARShip-1.0-1_r02432_c02432.tiff') as dataset:
    # 获取仿射变换矩阵
    transform = dataset.transform
    print("Transform:", transform)

    # 获取地理坐标参考系统（CRS）
    crs = dataset.crs
    print("CRS:", crs)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
Transform: | 3.00, 0.00, 398197.03|
| 0.00,-3.00, 117781.13|
| 0.00, 0.00, 1.00|

CRS: PROJCS["Projection Name = UTMUnits = Meter",

GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257164354493]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],

PROJECTION["Transverse_Mercator"],

PARAMETER["latitude_of_origin",30.1750794778181],

PARAMETER["central_meridian",114.404640526027],

PARAMETER["scale_factor",0.9996],

PARAMETER["false_easting",400000],

PARAMETER["false_northing",100000],

UNIT["metre",1,AUTHORITY["EPSG","9001"]],

AXIS["Easting",EAST],

AXIS["Northing",NORTH]]

第一个仿射矩阵是数据集自带的，但是第二仿射矩阵需要使用 pyporj 库种的 Transformer 进行计算，如下面代码所示，获得的仿射矩阵就是从数据集中的 crs 坐标系到“EPSG:4326”的映射。EPSG:4326 是经纬度的坐标系标准。

1
2
3
4
5
6
7
8
9
10
11
12
13
from pyproj import Transformer
crs = dataset.crs  #获取数据集的CRS
geo_transformer = Transformer.from_crs(crs, "EPSG:4326", always_xy=True) #获得将CRS转化为标准经纬度的坐标
def pixel_to_lonlat(x_pixel, y_pixel, transform_affine, geo_transformer):
    """Convert pixel coordinates to lon/lat using raster transform/CRS."""
    if transform_affine is None:
        raise ValueError("TIFF file does not contain an affine transform matrix.")
    x_geo, y_geo = transform_affine * (x_pixel, y_pixel)
    if geo_transformer is None:
        # CRS missing; return projected coordinates as best effort.
        return x_geo, y_geo
    lon, lat = geo_transformer.transform(x_geo,y_geo)
    return lon, lat

验证结果

使用公开的数据集 AIR-SARShip-1.0 进行验证，原始的图像如下所示：

image-20251124102828537

将获得的中心经纬度在 google 地球中进行搜索，结果如下，验证了经纬度的真实性。

image-20251124102816453

SAR 图像裁剪

如上面输出的仿射矩阵：

1
2
3
Transform: | 3.00, 0.00, 398197.03|
		  | 0.00,-3.00, 117781.13|
		  | 0.00, 0.00, 1.00           |

对于一个坐标点（x，y），仿射转化的结果是 3 _ x +398197.03, -3 _ y+ 119981.13,其中的 398197.03 和 117781.13 相当于是偏移量，是图像左上角的 CRS 坐标。3 和-3 相当于是缩放变量，是一个像素点单位对应的 CRS 坐标的单位。

所以如果要对图像进行裁剪的话，需要将偏移量进行修改，才能保证子图的仿射矩阵正确。

1
2
window = Window(col, row, w, h)
new_transform_matrix = rasterio.windows.transform(window, dataset.transform)

可以通过上面的方法获得裁剪之后的仿射矩阵，其中 window 是裁剪之后的图像在原图像中的框图。