首页 > 编程语言 >用Python计算栅格数据的真实面积

用Python计算栅格数据的真实面积

时间:2024-11-11 11:56:38浏览次数:4  
标签:真实 Python chunk 栅格数据 coords np input tif size

用Python计算栅格数据的真实面积

在地理空间分析中,栅格数据的像素值通常代表某种属性,比如土地利用比例、植被覆盖率等。这些数据往往基于经纬度网格表示的比例值,而为了更直观地理解这些数据的空间意义,我们需要将这些比例值转化为实际面积(如平方米或公顷)。

对于高分辨率的大尺寸栅格文件(GeoTIFF),一次性加载整个文件会耗费大量内存,甚至导致处理失败。为了解决这个问题,我们可以通过分块处理 的方式,逐步计算像素的真实面积,确保即使是大文件也能顺利处理。

一. 应用场景

该方法尤其适用于以下场景:

1. 比例值栅格数据:

当栅格像素值表示某个区域的比例(例如森林覆盖率)时,通过面积计算,可以得出真实的覆盖面积。

2. 土地利用分析:

计算农田、草地、建筑用地等各类土地利用类型的真实面积。

3. 遥感影像处理:

在遥感影像中将像素值的属性映射到真实的地理空间。

二.代码实现

import sys
import os
import rasterio
import numpy as np
from glob import glob

def haversine(coord1, coord2):
    """
    Calculate the great circle distance in kilometers between two points
    on the earth (specified in decimal degrees, with coordinates in the order of longitude, latitude).

    Arguments:
    coord1 -- tuple containing the longitude and latitude of the first point (lon1, lat1)
    coord2 -- tuple containing the longitude and latitude of the second point (lon2, lat2)

    Returns:
    distance in kilometers.
    """
    # Extract longitude and latitude, then convert from decimal degrees to radians
    lon1, lat1 = np.radians(coord1)
    lon2, lat2 = np.radians(coord2)

    # Haversine formula
    dlat = abs(lat2 - lat1)
    dlon = abs(lon2 - lon1)

    a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    c = 2 * np.arcsin(np.sqrt(a))
    r = 6371  # Radius of earth in kilometers. Use 3956 for miles
    return c * r


# Read geotiff file in data folder
# Calculate area for a single GeoTIFF file
def calculate_area_for_tif(input_tif, chunk_size=1024):
    # Check if input file exists
    if not os.path.isfile(input_tif):
        raise FileNotFoundError(f"Input file '{input_tif}' does not exist.")

    # Construct output file path
    output_path = input_tif.split('.')[0] + '_area.tif'

    # Skip if the area raster already exists
    if os.path.exists(output_path):
        print(f"Area of {input_tif} already calculated. Skipping...")
        return

    # Open the source raster
    with rasterio.open(input_tif) as src:
        # Check if the raster is in a geographic coordinate system
        if not src.crs.is_geographic:
            raise ValueError(f"The raster {input_tif} is not in a geographic coordinate system!")

        # Get raster metadata and size
        meta = src.meta.copy()
        width, height = src.width, src.height
        transform = src.transform

        # Update metadata for output
        meta.update(compress='lzw', dtype=rasterio.float32, count=1)

        # Calculate the total number of chunks
        total_chunks = ((height + chunk_size - 1) // chunk_size) * ((width + chunk_size - 1) // chunk_size)
        chunk_count = 0  # Track processed chunks

        # Create the output raster file
        with rasterio.open(output_path, 'w', **meta) as dst:
            print(f"Processing {input_tif} in chunks...")

            for row_start in range(0, height, chunk_size):
                row_end = min(row_start + chunk_size, height)

                for col_start in range(0, width, chunk_size):
                    col_end = min(col_start + chunk_size, width)

                    # Read a chunk of the source raster
                    window = rasterio.windows.Window(
                        col_start, row_start, col_end - col_start, row_end - row_start
                    )
                    chunk_transform = src.window_transform(window)

                    rows, cols = np.meshgrid(
                        np.arange(row_start, row_end),
                        np.arange(col_start, col_end),
                        indexing='ij'
                    )

                    # Apply geotransform to get geographical coordinates
                    x_coords, y_coords = chunk_transform * (cols, rows)
                    x_coords_right, y_coords_right = chunk_transform * (cols + 1, rows)
                    x_coords_bottom, y_coords_bottom = chunk_transform * (cols, rows + 1)

                    xy_coords = np.stack((x_coords, y_coords), axis=0)
                    xy_coords_right = np.stack((x_coords_right, y_coords_right), axis=0)
                    xy_coords_bottom = np.stack((x_coords_bottom, y_coords_bottom), axis=0)

                    # Calculate the area for this chunk
                    length_right = haversine(xy_coords, xy_coords_right)
                    length_bottom = haversine(xy_coords, xy_coords_bottom)
                    area_chunk = length_right * length_bottom

                    # Write the chunk to the output raster
                    dst.write(area_chunk.astype(np.float32), 1, window=window)

                    # Update progress
                    chunk_count += 1
                    progress = (chunk_count / total_chunks) * 100
                    print(f"Progress: {progress:.2f}% ({chunk_count}/{total_chunks} chunks)")

            print(f"\nArea of {input_tif} calculated and saved to {output_path}")


if __name__ == '__main__':
    # Specify the input file path
    input_tif = 'path/to/your/input_file.tif'

    # Run the function
    calculate_area_for_tif(input_tif, chunk_size=1024)
    print('\nArea calculation complete!')

! 本代码修改自迪肯大学王金柱博士的代码

标签:真实,Python,chunk,栅格数据,coords,np,input,tif,size
From: https://www.cnblogs.com/skypanxh/p/18539433

相关文章

  • python中常见的8种数据结构之一字典及其使用方法
    字典(Dictionary)是Python中常见的数据结构之一,用于存储一组配对的键(key)和值(value)。字典是可变的、无序的,并且键必须是唯一的。创建字典的方法有两种:使用花括号{}或使用内置的dict()函数。下面是一些常见的字典操作和方法:1.创建字典:my_dict={'key1':'value1','key2'......
  • 使用wxpython开发跨平台桌面应用,动态工具的创建处理
    在我们开发系统的时候,往往需要一个很容易理解功能的工具栏,工具栏是一个系统的快速入口,美观易用的工具栏是可以给系统程序增色不少的,本篇随笔介绍在使用wxpython开发跨平台桌面应用,工具栏的动态展现效果,以及多级工具栏显示等的创建处理过程。1、wxpython工具栏介绍在wxPython中......
  • python如何读取json文件
    在Python中读取JSON文件通常使用json模块,这是Python标准库的一部分,不需要额外安装。以下是读取JSON文件的基本步骤:打开JSON文件。使用json.load()函数将文件内容解析为Python对象(通常是字典或列表)。关闭文件。下面是一个具体的例子:pythonimportjson打开JSON文件withop......
  • Python 列表:数据处理的强大工具
    文章目录一、Python列表的基本概念二、Python列表的特性三、Python列表的操作方法四、Python列表在数据处理中的应用Python列表:数据处理的强大工具而在Python的众多数据结构中,列表(List)无疑是使用最为广泛的一种。一、Python列表的基本概念Python列表是一种有......
  • 京东商品详情,Python爬虫的“闪电战”
    在这个数字化的时代,我们每天都在和数据打交道,尤其是电商数据。想象一下,你是一名侦探,需要快速获取京东上某个商品的详细信息,但是没有超能力,怎么办?别担心,Python爬虫来帮忙!今天,我们就来一场幽默的“闪电战”,用Python快速获取京东商品详情。为什么选择Python做“武器”?选择Pytho......
  • Python爬虫快速获取JD商品详情:代码示例与技巧解析
    在当今这个信息爆炸的时代,数据成为了一种宝贵的资源。对于电商行业来说,获取商品详情信息是进行市场分析、价格比较、库存管理等重要环节的基础。本文将通过一个Python爬虫示例,展示如何快速获取(JD)商品的详情信息。为什么选择Python进行爬虫开发?Python作为一种高级编程语言,以......
  • 基于Python的大模型学习手册(基础级)
    前言大模型(全称为大语言模型,英文名称:LargeLanguageModel),这个2023年刷爆了互联网圈的“现象级明星”,几乎以前所未有的姿态,席卷了各行各业,世人一时为之惊叹。同时,也开辟了各大厂商投入AI研发的新赛道。前排提示,文末有大模型AGI-CSDN独家资料包哦!乘着这波“西风”,国内大......
  • python 制作智慧课堂点名系统
    #Python制作智慧课堂点名系统##一、项目背景-智慧课堂需求-点名系统的重要性##二、技术选型-Python语言介绍-适合的Python库(如tkinter,pandas等)##三、系统设计###3.1功能需求-学生名单管理-随机点名-点名记录保存与查询###3.2数据库设计-数据库......
  • 毕业设计:python考研院校推荐系统 混合推荐 协同过滤推荐算法 爬虫 可视化 Django框架(
    毕业设计:python考研院校推荐系统混合推荐协同过滤推荐算法爬虫可视化Django框架(源码+文档)✅1、项目介绍技术栈:Python语言MySQL数据库Django框架协同过滤推荐算法requests网络爬虫pyecharts数据可视化html页面、爬取院校信息:https://yz.chsi.com.cn/sch/(研招网......
  • 鸿蒙Next设备认证机制:Device Certificate Kit的真实性证明应用
    本文旨在深入探讨华为鸿蒙HarmonyOSNext系统(截止目前API12)的技术细节,基于实际开发实践进行总结。主要作为技术分享与交流载体,难免错漏,欢迎各位同仁提出宝贵意见和问题,以便共同进步。本文为原创内容,任何形式的转载必须注明出处及原作者。在当今数字化的浪潮中,设备的安全......