python逐像素获取栅格经纬度分别保存在两个矩阵中(代码收藏)

# -*- coding: utf-8 -*-
"""
任何问题联系邮箱：
chinesevoice@163.com
"""
from osgeo import gdal
import numpy as np
from itertools import product
from functools import partial
# import multiprocessing as mp
import time
gdal.AllRegister
dataset = gdal.Open('218-2020-07-01-2020-10-01.tif',gdal.GA_ReadOnly)
trans = dataset.GetGeoTransform()
proj = dataset.GetProjection()
rows = dataset.RasterYSize
cols = dataset.RasterXSize
print('坐标六参数:',trans)
print('投影参数:',proj)
print('行列数：',rows,cols)
#==========第一种方法===========
#行列转经纬度
def rowcol2lonlat(extend,xsize,ysize):
    #xsize表示列
    lon = extend[0] + xsize * extend[1] + ysize * extend[2] + extend[1]/2
    lat = extend[3] + xsize * extend[4] + ysize * extend[5] + extend[5]/2
    return lon,lat
#分别获取两个矩阵
def get_lon_lat_array(row,col,tran):
    lon_array = np.repeat(np.nan, row*col).reshape(row,col)
    lat_array = np.repeat(np.nan, row*col).reshape(row,col)
    for r in range(0,row):
        for c in range(0,col):
            lon_array[r,c],lat_array[r,c] = rowcol2lonlat(tran,c,r)
    return lon_array,lat_array
ts = time.time()
lon_array,lat_array = get_lon_lat_array(rows,cols,trans)
# print(lon_array.shape,lat_array.shape)
# print(lon_array.view())
te = time.time()
print('第一种方法用时：',te-ts)
#==========第二种方法===========
def get_lon_lat_array_2(row,col,tran):
    rowlst = np.arange(0, row,step = 1)
    collst = np.arange(0, col,step = 1)
    #rc_comb = product(rowlst,collst)
    rc2lon = lambda rowcol,tran:tran[0] + rowcol[1] * tran[1] + rowcol[0] * tran[2] + tran[1]/2
    rc2lat = lambda rowcol,tran:tran[3] + rowcol[1] * tran[4] + rowcol[0] * tran[5] + tran[5]/2
    # lon_array = np.array(list(map(partial(rc2lon,tran = trans),product(rowlst,collst)))).reshape(row,col)
    # lat_array = np.array(list(map(partial(rc2lat,tran = trans),product(rowlst,collst)))).reshape(row,col)
    # lon_array = np.fromiter(map(partial(rc2lon,tran = trans),product(rowlst,collst)),dtype=np.float64).reshape(row,col)
    # lat_array = np.fromiter(map(partial(rc2lat,tran = trans),product(rowlst,collst)),dtype=np.float64).reshape(row,col)
    lon_array = [*map(partial(rc2lon,tran = trans),product(rowlst,collst))]
    lat_array = [*map(partial(rc2lat,tran = trans),product(rowlst,collst))]
    return lon_array,lat_array
ts2 = time.time()
lon_array2,lat_array2 = get_lon_lat_array_2(rows,cols,trans)
# print(lon_array2.shape,lat_array2.shape)
# print(lon_array2.view())
te2 = time.time()
print('第二种方法用时：',te2-ts2)