#%% 读取文本文件到数据中
import numpy as np
s = open('no2_202009.asc').readlines() #>>>>>>>>>>>>>>>>>>>>>修改你的文件路径
newArray = np.zeros((1440,2880)) -9 #空的数组，并设置noData为-9
# 145: 每145行是一组，线标注了纬度值,然后是该纬度下所有的值
# 读取每一个纬度下的数据，跳过前 4+1 行说明
for i in range(4,len(s),145):
# print(s[i]) #输出纬度值
sTmp = [] #列表置空
for l in s[i+1:i+145]: #同一纬度下各个经度所对应的值
l=l.replace('\n','') #去除末尾的转义字符
sTmp.append(l) #拼到大的列表中
sTmp = ''.join(sTmp) #拼接为一个大的字符串
i_each = (i-4)//145
for j in range(0,len(sTmp),4): #每4个字符是一个数
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
#%% 将数组存储为 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='')
以下内容为 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()
因为第一行的数据是-90°的，而我们希望矩阵的左下角才是，所以需要调换过来，但是，本来应该更容易解决的问题可能由于字符编码的问题总是出错，干脆直接用matlab写个简单的行调换。
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
之前还写了一个中国区域的
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: #东经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_each+1,jj,int(s[i][j:j+4]))
print(i)
f.close()
