#%% 讀取文本文件到數(shù)據(jù)中
import numpy as np

s = open('no2_202009.asc').readlines() #>>>>>>>>>>>>>>>>>>>>>修改你的文件路徑

newArray = np.zeros((1440,2880)) -9 #空的數(shù)組，并設(shè)置noData為-9

# 145: 每145行是一組，線標(biāo)注了緯度值,然后是該緯度下所有的值
# 讀取每一個(gè)緯度下的數(shù)據(jù)，跳過(guò)前 4+1 行說(shuō)明
for i in range(4,len(s),145): 
    # print(s[i]) #輸出緯度值
    sTmp = [] #列表置空
    for l in s[i+1:i+145]: #同一緯度下各個(gè)經(jīng)度所對(duì)應(yīng)的值
        l=l.replace('\n','') #去除末尾的轉(zhuǎn)義字符
        sTmp.append(l) #拼到大的列表中
    sTmp = ''.join(sTmp) #拼接為一個(gè)大的字符串
    
    i_each = (i-4)//145
    for j in range(0,len(sTmp),4): #每4個(gè)字符是一個(gè)數(shù)
        j_each = j//4
        
        if int(sTmp[j:j+4])>0:
            newArray[i_each][j_each] = int(sTmp[j:j+4])
        else:
            newArray[i_each][j_each] = -9


#%% 將數(shù)組存儲(chǔ)為 nc 或 tif
import xarray as xr
newArray = np.where(newArray==-9,np.nan,newArray)
latArray = np.array(range(1440))*0.125-89.9375
longArray = np.array(range(2880))*0.125-179.9375

ds = xr.Dataset( {'tropNO2': (("latitude", "longitude"), newArray )},
                coords= { 
                    'latitude': latArray, 
                    'longitude': longArray
                        }
               )
ds['tropNO2'].attrs = {'units':'10^13 molec. cm-2'}

ds.to_netcdf('no2_202009.nc')


#%% 也可存儲(chǔ)為 tif 格式
import rioxarray
ds = ds.rio.write_crs("epsg:4326")
ds.rio.to_raster("no2_202009.tif")


#%% 顯示圖像
# import matplotlib.pyplot as plt
# plt.imshow(newArray[::-1]) #簡(jiǎn)單繪制出影像


#%% 存儲(chǔ)為 ArcGIS 或 QGIS 可以打開(kāi)的格式
# headerArcGIS = '''ncols 2880
# nrows 1440
# xllcenter -179.9375
# yllcenter -89.9375
# cellsize 0.125
# NODATA_value -9'''
# np.savetxt("no2_202009_new.asc", newArray[::-1], fmt="%d", delimiter=" ",header=headerArcGIS,comments='') 

以下內(nèi)容為 2018 年寫(xiě)的

import numpy as np
s = open ('C:\\Users\\tropomi\\no2_201810.asc').readlines()
print(len(s))
a = np.zeros((1440,2880))
#f = open('10.txt','w+')
for i in range(4,len(s)):
    ii = (i-4)//145*0.125-89.9375
    i_each = (i-4)%145
    for j in range(0,80,4):
        if i_each>0 and int(s[i][j:j+4])>0:
            j_each = j//4
            jj = ((i_each-1)*20+j_each)*0.125-179.9375
            #print(ii,jj,s[i][j:j+4],file=f)
            a[(i-4)//145][(i_each-1)*20+j_each] = int(s[i][j:j+4])
            #a[1339-((i-4)//145)][(i_each-1)*20+j_each] = int(s[i][j:j+4])#可能字符編碼的問(wèn)題
    print(i)
np.savetxt("10.txt", a, fmt="%d", delimiter=" ")
#f.close()

因?yàn)榈谝恍械臄?shù)據(jù)是-90°的，而我們希望矩陣的左下角才是，所以需要調(diào)換過(guò)來(lái)，但是，本來(lái)應(yīng)該更容易解決的問(wèn)題可能由于字符編碼的問(wèn)題總是出錯(cuò)，干脆直接用matlab寫(xiě)個(gè)簡(jiǎn)單的行調(diào)換。

clear;
clc;
filename = '.\02.txt';
a=textread(filename);
b=a;
for i=1:1440
    b(1441-i,:)= a(i,:);
end

dlmwrite('02new.txt', b, 'delimiter', ' ','precision', 5,'newline', 'pc')

最后在txt的頭部添上信息，可參照ArcMap的說(shuō)明

ncols 2880
nrows 1440
xllcenter -179.9375
yllcenter -89.9375
cellsize 0.125
NODATA_value 0

之前還寫(xiě)了一個(gè)中國(guó)區(qū)域的

s = open('/Users/heqin/Downloads/tropomi/no2_201809.asc').readlines()
f = open('9.txt','w+')
for i in range(107885,167044): #北緯3度～54度
    ii = (i-4)//145*0.125-89.9375
    i_each = (i-4)%145 #該i在該組的實(shí)際第幾行
    for j in range(0,80,4):
        if 102 <= i_each <= 127 and int(s[i][j:j+4]) >0: #東經(jīng)73度～136度
            j_each = j//4 #該j在該組的實(shí)際第幾列
            jj = ((i_each-1)*20+j_each)*0.125-179.9375
            print(ii,jj,int(s[i][j:j+4])/100,file = f) #單位換成10^15
            #print("實(shí)際行和緯度",i+1,ii,"該行的第X個(gè)的經(jīng)度",j_each+1,jj,int(s[i][j:j+4]))
    print(i)
f.close()