Python nc文件批量转tif再转投影

import warnings

import nc
import netCDF4
warnings.filterwarnings('ignore')
warnings.filterwarnings('ignore', category=DeprecationWarning)
import netCDF4
import pandas as pd
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt
import math
import h5py
from pyhdf.SD import *
from osgeo import osr
import numpy as np

###
def array2raster(newRasterfn, rasterOrigin, xsize, ysize, array):
    """
     newRasterfn: 输出tif路径
     rasterOrigin: 原始栅格数据Extent
     xsize: x方向像元大小
     ysize: y方向像元大小
     array: 计算后的栅格数据
    """
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    originX = rasterOrigin[0]  # 起始像元经度
    originY = rasterOrigin[1]  # 起始像元纬度
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    # 括号中两个0表示起始像元的行列号从(0,0)开始
    outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
    # 获取数据集第一个波段,是从1开始,不是从0开始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    # 代码4326表示WGS84坐标
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


###
from os import listdir
import re

inputf = 'E:/FSC_AVHRR/train_nc/'
output_TIF = 'E:/FSC_AVHRR/train_tif_4326/'
output_TIF_RE = 'E:/FSC_AVHRR/train_tif_6933/'
files = listdir(inputf)

pattern = '.nc$'
timer = 0

for i in files:
    temp_input = inputf + i

    temp_out_name = re.split(pattern, i)[0] + '.tif'
    temp_out_T_name = re.split(pattern, i)[0] + '_T.tif'
    temp_output = output_TIF + temp_out_name
    temp_output_T = output_TIF_RE + temp_out_T_name

    nc_data = netCDF4.Dataset(temp_input)
    FSC = nc_data.variables['scfg'][:]
    # plt.imshow(C)
    FSC_data = np.flipud(FSC[0, :, :])

    xsize = 0.05
    ysize = -0.05
    # 对于全球尺度rasterOrigin 为[-180(xmin),-90(ymin)]
    array2raster(temp_output, [-180, 90], xsize, ysize, FSC_data)
    option4 = gdal.WarpOptions(format='GTiff', srcSRS='EPSG:4326', dstSRS='EPSG:6933')
    FSC_RE = gdal.Warp(temp_output_T, temp_output, options=option4)
    timer = timer + 1
    print(timer)

 

posted @ 2022-04-15 16:22  揪你小辫子  阅读(479)  评论(0编辑  收藏  举报