How to create a netCDF file with python netCDF4 from tif?

Question

I would like to convert a folder with tif image to netcdf

import numpy as np
import datetime as dt
import os
import gdal
import netCDF4
import re
import xarray as xr

#PRIMERA IMAGEN 
ds = gdal.Open("first_image")

a = ds.ReadAsArray()
nlat, nlon = np.shape(a)

b = ds.GetGeoTransform()  # bbox, interval
lon = np.arange(nlon)*b[1]+b[0]
lat = np.arange(nlat)*b[5]+b[3]

#FECHA DE LA PRIMERA IMAGEN
basedate = dt.datetime(2020,10,16,9,0,0)

# Create NetCDF file
nco = netCDF4.Dataset("File.nc",
                      'w', format='NETCDF4')

nco.description = 'Example simulation data'

# Create dimensions, variables and attributes:
nco.createDimension('lon', nlon)
nco.createDimension('lat', nlat)
nco.createDimension('time', None)

timeo = nco.createVariable('time', 'f4', ('time',))
timeo.units = 'hours since 2012-09-28 00:00:00'
timeo.standard_name = 'time'

lono = nco.createVariable('lon', 'f4', ('lon',))
lono.standard_name = 'longitude'

lato = nco.createVariable('lat', 'f4', ('lat',))
lato.standard_name = 'latitude'

# Create short integer variable for temperature data, with chunking
tmno = nco.createVariable('tmn', 'f4',  ('time', 'lat', 'lon'),
                          zlib=True, 
                          fill_value=-9999)
tmno.units = 'degC'
tmno.scale_factor = 1#0.01
tmno.add_offset = 0.00
tmno.long_name = 'minimum monthly temperature'
tmno.standard_name = 'air_temperature'
tmno.set_auto_maskandscale(False)

nco.Conventions = 'CF-1.6'

# Write lon,lat
lono[:] = lon
lato[:] = lat


itime = 0

# Step through data, writing time and data to NetCDF
for root, dirs, files in os.walk(r'folder_images'):
    dirs.sort()
    files.sort()
    for f in files:
        if f[-3:]=='tif':
            print(f)
            # read the time values by parsing the filename
            year = 2012
            mon = 9
            day=28
            hora=int(f[0:2])
            date = dt.datetime(year, mon, day, hora, 0, 0)
            print(date)
            dtime = (date-basedate).total_seconds()/3600
            timeo[itime] = dtime
            # min temp
            tmn_path = os.path.join(root, f)
            print(tmn_path)
            tmn = gdal.Open(tmn_path)
            a = tmn.ReadAsArray()  # data
            tmno[itime, :, :] = a
            itime = itime+1
            tmn=None
nco.close()
ds =None

When I add the netcdf file to Qgis the coordinates are wrong:

Extent 0.5000000000000000,0.5000000000000000 : 804.5000000000000000,559.5000000000000000

The raster extent are: Extent 576172.0000000000000000,4143175.8648832864128053 : 616453.8192302760435268,4171198.0000000000000000

what am I doing wrong?

Images: https://drive.google.com/file/d/13TkWx3UgSyBt48nm0UMC6GJRO_Nwoch6/view?usp=sharing

Answer 1

等待大神答复