CAMS再分析数据逐日提取

前注：
在处理CAMS再分析资料时我下载了一年的数据，官网直接把一年的数据直接打包成了一个nc给我，对于从没见到过这种世面的我，当场就傻了，然后就想着给他提取出来，其实原理很简单，把这个三维的数组给一层一层剥离出来，伪代码如下：

循环三维数组的层：
  读取出当前层的数据
  读出对应层对应的时间
  读取经纬度信息以及行列信息
  按照上面的信息给他写成常用于处理的tif文件
结束循环

但是事情总是不会让你那么顺心，马上就遇到了第一个问题：时间怎么转换，这个时间是这样的：

这个是按照小时然后从1900年1月1日的00时开始记录的，那也没问题，我们转就行了，我想到闰年和平年的区别，不想再多写那几行，所以就偷懒，打算直接用python自带的time和datetime包完成 时间戳的转换。结果发现，python的时间戳转换最早时间只支持到1970年的1月1号。可恶啊！！
找了好久，那相当于python中的1970年时0，那我们直接把之前的总小时算出来然后再结合时间戳不就行了，具体思路参考了博客：
python时间戳转换（1970前）
前面加上netcdf包读取文件的方式，将这个文件先读取出来，然后就以时间轴为他的层数吧，每层代表的一个时间，按照时间先后设置的第三维：
然后就是计算我们1970年前的总共时间：这里注意：1900年是平年，而不是闰年，一定要搞清楚，不然就怎么都不能给他把时间换对

time_total=0.0
for i in range(1900,1970):
    year=int(i)#得到当前年份
    if year % 4 ==0:#如当前年份能被四整除，那这一年就是闰年，有366天
        year_sec=366*24
    else:#否则是平年，只有365天
        year_sec=365*24
    time_total+=year_sec
time_total=time_total-24#1900年能被100整除，所以也是平年（百年不润）

下面记录他的经纬度信息：其实能够直接知道它的经纬度的分辨率就是0.75°，所以我也不做过多解释了：

len_lon=int(360/0.75)#经度维
len_lat=int(180.75/0.75)#纬度维
latmin=-90
latmax=90
lonmin=-180
lonmax=180#其实这个的经纬度最大最小不是这样的，这儿就不用改了，我在前面改了（180.75那里，不然他的格点数和图像对不上）
resolution=0.75

后边就是直接出图，这里面其实主要用到的就是netcdf4、datetime、time、osgeo的gdal和osr以及numpy包

import netCDF4 as nc
import datetime
import time as tt
from osgeo import gdal
from osgeo import osr
import numpy as np
filename="F:\\CAMS_tif\\2020.nc"
f=nc.Dataset(filename)#如果对netcdf有不懂的可以搜其他的博客，我前面的fnl再分析资料里面也提到了一点点
tcch4=f.variables['tcch4']
all_vars_info=f.variables.items()
# print(all_vars_info)
time=f.variables['time']
nc_time=time[:,]
#time begins: 1900 01 01 00:00:00.0
time_total=0.0
for i in range(1900,1970):
    year=int(i)
    if year % 4 ==0:
        year_sec=366*24
    else:
        year_sec=365*24
    time_total+=year_sec
time_total=time_total-24
# print(time_total)
# print(datetime.datetime(1970,1,1)+datetime.timedelta(hours=-(time_total)))
# 下面用nc的时间先减去1970前的时间戳
len_lon=int(360/0.75)
len_lat=int(180.75/0.75)
latmin=-90
latmax=90
lonmin=-180
lonmax=180
resolution=0.75
for i in range(len(nc_time)):
    time_i=nc_time[i]
    time_delta=int(time_i)-time_total
    time_now=datetime.datetime(1970,1,1,0,0)+datetime.timedelta(hours=(time_delta))#这儿主要就是时间戳的转换，前面计算了70年前的delta时间，这儿直接加起来就行啦！
    time_now=str(time_now)
    time_year=time_now[0:4]
    time_month=time_now[5:7]
    time_day=time_now[8:10]
    time_hour=time_now[11:13]
    file_date=time_year+time_month+time_day+time_hour
    # print(file_date)
    driver=gdal.GetDriverByName('GTiff')
    out_tif_name=r'F:/CAMS_tif/CH4/{}{}{}'.format('CH4_',file_date,'.tiff')
    out_tif=driver.Create(out_tif_name,len_lon,len_lat,1,gdal.GDT_Float32)
    geotransfor=(lonmin,resolution,0,latmax,0,-resolution)
    out_tif.SetGeoTransform(geotransfor)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)#WGS84
    out_tif.SetProjection(srs.ExportToWkt())
    variable_out=tcch4[i]
    data_trans=np.asarray(variable_out)
    out_tif.GetRasterBand(1).WriteArray(data_trans)
    out_tif.FlushCache()
    out_tif=None
    print('The file:'+file_date+' has been done!')
print('*****************************************************')
print('All process has been done!')
print('*****************************************************')