锐多宝
锐多宝的朋友圈
锐多宝

//// 批量对齐栅格数据代码

最近经常遇到把一个栅格数据对齐到另外一个栅格数据。

经常会遇到同一个区域的几个栅格数据的分辨率、行列号和投影都不一致。

那可以用gdal和rasterio来实现。代码如下:

import rasterio
from rasterio.warp import reproject, Resampling
import os

def align_raster(input_path, reference_path, output_path):
    """
    将输入栅格对齐到参考栅格
    """
    try:
        // 读取参考栅格获取目标参数
        with rasterio.open(reference_path) as src2:
            dst_crs = src2.crs
            dst_transform = src2.transform
            dst_width = src2.width
            dst_height = src2.height

        // 对输入栅格进行重投影和重采样
        with rasterio.open(input_path) as src1:
            // 准备输出元数据
            meta = src1.meta.copy()
            meta.update({
                'crs': dst_crs,
                'transform': dst_transform,
                'width': dst_width,
                'height': dst_height
            })
            
            // 创建输出文件并进行重投影
            with rasterio.open(output_path, 'w', **meta) as dst1:
                for i in range(1, src1.count + 1):
                    reproject(
                        source=rasterio.band(src1, i),
                        destination=rasterio.band(dst1, i),
                        src_transform=src1.transform,
                        src_crs=src1.crs,
                        dst_transform=dst_transform,
                        dst_crs=dst_crs,
                        resampling=Resampling.bilinear
                    )
        print(f"已保存对齐后的栅格: {output_path}")
        return True

    except Exception as e:
        print(f"处理文件 {input_path} 时发生错误: {str(e)}")
        return False

def main():
    // 定义输入输出路径
    input_dir = r'你的文件路径'
    output_dir = r'输出路径'
    reference_path = r'参考栅格的tif'

    // 确保输出目录存在
    os.makedirs(output_dir, exist_ok=True)

    // 获取所有tif文件
    tif_files = [f for f in os.listdir(input_dir) if f.endswith('.tif')]
    total_files = len(tif_files)
    
    print(f"共找到 {total_files} 个TIF文件需要处理")

    // 处理每个文件
    successful = 0
    failed = 0
    
    for i, tif_file in enumerate(tif_files, 1):
        input_path = os.path.join(input_dir, tif_file)
        output_path = os.path.join(output_dir, tif_file)
        
        print(f"\n处理第 {i}/{total_files} 个文件: {tif_file}")
        
        if align_raster(input_path, reference_path, output_path):
            successful += 1
        else:
            failed += 1

    // 打印处理结果统计
    print("\n处理完成:")
    print(f"成功处理: {successful} 个文件")
    print(f"处理失败: {failed} 个文件")
    print(f"总计文件: {total_files} 个")

    // 验证最后一个处理的文件
    if successful > 0:
        print("\n验证最后处理的文件:")
        with rasterio.open(output_path) as aligned1, rasterio.open(reference_path) as src2:
            print(f"参考栅格大小: {src2.width} x {src2.height}")
            print(f"对齐后栅格大小: {aligned1.width} x {aligned1.height}")
            print(f"参考栅格变换矩阵: {src2.transform}")
            print(f"对齐后变换矩阵: {aligned1.transform}")
            
            if (aligned1.transform == src2.transform and 
                aligned1.width == src2.width and 
                aligned1.height == src2.height):
                print("栅格已成功对齐到参考栅格的网格和分辨率")
            else:
                print("对齐结果有误,请检查")

if __name__ == "__main__":
    main()

原理是:

  • 读取参考栅格,获取目标参数:

    • 坐标系统(CRS)
    • 变换矩阵(transform)
    • 影像宽度和高度
  • 对每个输入栅格文件:

    • 读取源文件
    • 更新元数据(metadata)
    • 使用双线性重采样方法进行重投影
    • 保存对齐后的文件

注:该代码由chatgpt o1一次性生成,生成后就可以直接运行,没有报错。

2024年11月15日
锐多宝

回头看landcover100的代码,不到一年,
感觉好陌生啊,这竟然是我自己敲的,卧槽。
真NB
修复一下这个网站。

2024年10月08日
锐多宝

SB很多,我也是其中一个。

2024年10月08日
锐多宝

一个通过照片识别地理位置的应用

引言

最近发现一个能根据照片进行地理位置判定的应用,在全球范围内能够非常准确地进行空间位置识别。我分3个尺度进行了测试,分别是城市街景(来源google和腾讯街景)、野外街景和我自己拍摄的照片进行测试。

街景图片(城市)测试

1.孟买街景测试

实际:北纬19°,经度72°

预测:北纬19°,经度72°

解释:这张照片拍摄于印度孟买的火柴厂巷。这张照片是一条两边都有几栋建筑的街道。有人走在街上,有几棵树。这些建筑大多是混凝土建造的,有很多窗户。这条街很窄,车流量很大。

2.伦敦街景测试

实际:北纬51°,经度0°

预测:北纬51°, 经度0°

3.迪拜街景测试

实际:北纬55°, 东经25°

预测:北纬55°, 东经25°

解释:这张照片是在阿拉伯联合酋长国迪拜拍摄的。这可以由背景中建筑的建筑来决定,这是城市的特征。此外,路标是阿拉伯联合酋长国的官方语言阿拉伯语。

4.悉尼街景测试

实际:南纬33°, 东经151°

预测:南纬33°, 东经151°

解释:这张照片是在澳大利亚新南威尔士州的阿什菲尔德拍摄的。这是由以下线索决定的:
-这些建筑的结构与阿什菲尔德的一致。
-植被与阿什菲尔德一致。
-标识与阿什菲尔德的标识一致。

5.旧金山街景测试

实际:北纬37°, 西经122°

预测:北纬37°, 西经122°

解释:这张照片是在旧金山的一座高楼上拍摄的。这是城市市中心的景色,远处是海湾大桥。这张照片是在一个雾天拍摄的,这在旧金山很常见。

6.成都街景测试

实际:北纬30°, 东经104°

预测:北纬31°, 西经121°

解释:这张照片是在中国拍摄的。标志上的字是中文的。这些建筑的建筑也是中国式的。

7.南京街景测试

实际:北纬32°, 东经118°

预测:北纬31°, 西经121°

解释:这张照片是在中国拍摄的,因为标志上的文字是中文的。这张照片是在一个城市里拍摄的,因为那里有高楼大厦和很多交通。这张照片是在街上拍的,因为路上有车。

街景图片(乡村与野外)测试

1.日本乡村测试

实际:北纬35°, 东经138°

预测:北纬35°, 西经136°

解释:这张照片是在日本的一个农村地区拍摄的。路又窄又弯,树又高又密。唯一标志是远处的一所小房子。

2.克罗利亚荒野测试

实际:北纬46°, 东经18°

预测:北纬44°, 西经20°

解释:这张照片拍摄于塞尔维亚的一个农村地区。道路被田野包围,看不到任何建筑或其他结构。唯一可见的标志是指示限速为60公里/小时的路标。植被与塞尔维亚的植被一致,景观和建筑也是如此。

3.泰国乡村测试

实际:北纬14°, 东经101°

预测:北纬14°, 西经100°

解释:这张照片是在泰国的一个农村地区拍摄的。道路没有铺设,植被茂密。看不到任何建筑物或其他构筑物。

我自己拍摄的照片测试

实际:北纬45°, 东经126°

预测:北纬39°,东经116°

解释:这张照片是在中国拍摄的。标牌上的字是中文。背景中的树木是典型的中国风景。

实际:北纬18°, 东经109°

预测:北纬18°, 东经109°

解释:这张照片是在南中国海的一个码头上拍摄的。背景中的雕像是海南南山观音,一尊108米高的观音菩萨雕像。这座雕像位于海南省三亚市的南山岛上

总结

这个网站是我目前遇到过对于照片位置预测最准确的一个应用,不仅能在城区准确预测,也能在乡村比较准确。能根据背景文字、植被、道路情况进行综合判定。目前该应用完全免费使用且无需登录,并且提供API接口(网址:https://geospy.ai/)。

需要明确,我目前的测试方法并不科学,测试样本少且设计也不科学。但不可否认该应用能大致的判定全球范围内的照片拍摄位置

我没有从事这方面的研究,也想知道一下它的预测原理,以及为什么能达到如此的高精度。

参考:

https://geospy.ai/

2024年05月05日
锐多宝

使用gdal均匀筛选点矢量

作用:

通过计算各点之间的欧式距离,筛选出符合目标的、均匀发布在空间中的N个数据点。

效果示意图

运行环境

python 3.10
安装:tqdm、numpy和tqdm这三个库

完整代码

import numpy as np
from osgeo import ogr, osr
from tqdm import tqdm

# 代码作用:通过计算各点之间的欧式距离,筛选出符合目标的、均匀发布在空间中的N个数据点。



# 定义需要采样的个数
n_samples = 100
input_path = r"测试数据\村点.shp"
output_path = r"测试数据\samples.shp"

# 1. 读取原始点数据
driver = ogr.GetDriverByName('ESRI Shapefile')
inds = driver.Open(input_path, 0)
layer = inds.GetLayer()

# 2. 提取点坐标和属性
coords = []
attrs = []
for feature in layer:
    geom = feature.GetGeometryRef()
    coords.append((geom.GetX(), geom.GetY()))
    attrs.append([feature.GetField(i) for i in range(feature.GetFieldCount())])
coords = np.array(coords)
attrs = np.array(attrs)

# 3. 定义距离函数
def distance(p1, p2):
    return np.sqrt(np.sum((p1 - p2)**2))

# 4. 随机选择第一个点
idx = np.random.choice(coords.shape[0], 1)
samples = coords[idx]
sample_attrs = attrs[idx]

# 5. 选择空间均衡的采样点
for _ in tqdm(range(n_samples - 1)):
    dists = np.array([np.min(np.array([distance(p, s) for s in samples])) for p in coords])
    idx = np.argmax(dists)
    samples = np.append(samples, [coords[idx]], axis=0)
    sample_attrs = np.append(sample_attrs, [attrs[idx]], axis=0)
    coords = np.delete(coords, idx, axis=0)
    attrs = np.delete(attrs, idx, axis=0)

# 6. 将采样点转为gdal几何对象
out_samples = []
for sample in samples:
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(sample[0], sample[1])
    out_samples.append(point)

# 7. 创建新的矢量层并写入采样点
out_driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = out_driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer('samples', layer.GetSpatialRef(), ogr.wkbPoint)

# 添加属性字段
for i in range(len(layer.schema)):
    field_defn = layer.schema[i]
    out_layer.CreateField(field_defn)

# 写入采样点要素
for i, sample in enumerate(out_samples):
    feature = ogr.Feature(out_layer.GetLayerDefn())
    feature.SetGeometry(sample)
    for j, attr in enumerate(sample_attrs[i]):
        feature.SetField(j, attr)
    out_layer.CreateFeature(feature)

out_layer = None
out_ds = None
2024年04月22日
锐多宝 .