#%% 讀取文本文件到數(shù)據(jù)中
import numpy as np
s = open('no2_202009.asc').readlines() #>>>>>>>>>>>>>>>>>>>>>修改你的文件路徑
newArray = np.zeros((1440,2880)) -9 #空的數(shù)組，并設置noData為-9
# 145: 每145行是一組，線標注了緯度值,然后是該緯度下所有的值
# 讀取每一個緯度下的數(shù)據(jù)宪赶，跳過前 4+1 行說明
for i in range(4,len(s),145):
# print(s[i]) #輸出緯度值
sTmp = [] #列表置空
for l in s[i+1:i+145]: #同一緯度下各個經(jīng)度所對應的值
l=l.replace('\n','') #去除末尾的轉(zhuǎn)義字符
sTmp.append(l) #拼到大的列表中
sTmp = ''.join(sTmp) #拼接為一個大的字符串
i_each = (i-4)//145
for j in range(0,len(sTmp),4): #每4個字符是一個數(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ù)組存儲為 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')
#%% 也可存儲為 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]) #簡單繪制出影像
#%% 存儲為 ArcGIS 或 QGIS 可以打開的格式
# 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 年寫的
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])#可能字符編碼的問題
print(i)
np.savetxt("10.txt", a, fmt="%d", delimiter=" ")
#f.close()
因為第一行的數(shù)據(jù)是-90°的，而我們希望矩陣的左下角才是脯燃，所以需要調(diào)換過來搂妻，但是，本來應該更容易解決的問題可能由于字符編碼的問題總是出錯辕棚，干脆直接用matlab寫個簡單的行調(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的說明
ncols 2880
nrows 1440
xllcenter -179.9375
yllcenter -89.9375
cellsize 0.125
NODATA_value 0
之前還寫了一個中國區(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在該組的實際第幾行
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在該組的實際第幾列
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("實際行和緯度",i+1,ii,"該行的第X個的經(jīng)度",j_each+1,jj,int(s[i][j:j+4]))
print(i)
f.close()
