OCO-2的NC格式转tif

最新推荐文章于 2024-10-07 01:36:35 发布

Desley、

最新推荐文章于 2024-10-07 01:36:35 发布

阅读量351

点赞数 3

文章标签： python windows

本文链接：https://blog.csdn.net/m0_62076976/article/details/138028820

版权

本文介绍了如何将OCO-2卫星的NC格式数据，特别是条带型一维数组，通过先转换为CSV，再利用点转栅格方法转化为GIS兼容的TIFF格式，以便于进一步分析和可视化。

摘要由CSDN通过智能技术生成

最近在写毕业论文，处理OCO-2卫星数据，该数据格式是NC，与大多数NC数据不同，该数据是一维数组，而非格网数据。即经度、维度、数值，且该数据为条带型

直接说方法：先读取nc格式中的经纬度、值，写入csv文件中，随后读取csv，利用类似Arcgis中的点转栅格的方式转为tif文件；

（1）nc转csv

import numpy as np
import netCDF4 as nc
import csv

# nc文件的路径
filename = r"D:\oco2_LtCO2_210101_B11100Ar_230607162301s.nc4"
f = nc.Dataset(filename)  # 相当于打开文件，并读取到变量f中
var_lon = f['longitude'][:]
var_lat = f['latitude'][:]
data = f['xco2_x2019'][:]
data_arr = np.asarray(data)

# 创建CSV文件并写入数据
with open('f:/output.csv', 'w', newline='') as csvfile:
    writer = csv.writer(csvfile)
    writer.writerow(['Longitude', 'Latitude', 'Data'])
    for i in range(len(var_lon)):
        if 90 <= var_lon[i] <= 115 and 25<= var_lat[i] <= 35:
            writer.writerow([var_lon[i], var_lat[i], data_arr[i]])

（2）csv转tif

import numpy as np
import pandas as pd
from osgeo import gdal, osr
import os

os.environ['PROJ_LIB'] = r"E:\py39\Lib\site-packages\osgeo\data\proj"


def create_tiff(data, lon, lat, data_field, output_filename, pixel_size=0.05):
    # 获取经纬度范围
    min_lon, max_lon = lon.min(), lon.max()
    min_lat, max_lat = lat.min(), lat.max()

    # 计算栅格的行列数
    num_cols = int((max_lon - min_lon) / pixel_size)
    num_rows = int((max_lat - min_lat) / pixel_size)

    # 创建一个空的栅格数组
    raster_data = np.zeros((num_rows, num_cols))

    # 将点数据分配到栅格单元中
    for i in range(len(data)):
        row = int((max_lat - lat[i]) / pixel_size) - 1
        col = int((lon[i] - min_lon) / pixel_size) - 1
        raster_data[row, col] = data[i]
    raster_data[raster_data == 0] = np.nan
    # 将栅格数据写入TIFF文件
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(output_filename, num_cols, num_rows, 1, gdal.GDT_Float32)
    dataset.SetGeoTransform((min_lon, pixel_size, 0, max_lat, 0, -pixel_size))

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # WGS84 EPSG code
    dataset.SetProjection(srs.ExportToWkt())

    band = dataset.GetRasterBand(1)
    band.WriteArray(raster_data)
    band.SetNoDataValue(-9999)
    band.FlushCache()

    dataset = None  # 关闭文件


# 读取CSV文件
csv_file = "F:\output2.csv"
data_df = pd.read_csv(csv_file)

# 提取所需的字段
lon = data_df['Longitude'].values
lat = data_df['Latitude'].values
data = data_df['Data'].values

# 设置输出文件名
output_filename = 'f:\output.tif'

# 调用函数创建TIFF文件
create_tiff(data, lon, lat, 'data值', output_filename)

print("TIFF文件已生成：", output_filename)

当然，也可以把两个代码合为一个，直接利用（2）中的函数处理（1）的数据即可