To describe my problem, I applied Polymer on my Sentinel-2 product over Valencia port, I have a NETCDF in output and I would like to recover each bands in GeoTIFF in lat/lon format.
I used gdal and python3 and the best way to have a georeferenced tiff, was the following commands:
Code: Select all
import sys, os, osgeo from osgeo import gdal import numpy as np from netCDF4 import Dataset from glob import glob as glb import subprocess filepath = sys.argv #netcdf file path filepathout = sys.argv #directory out f_in = Dataset(filepath, 'r') #to read netcdf lon = f_in.variables['longitude'][:] lat = f_in.variables['latitude'][:] Rw490 = f_in.variables['Rw490'][:] subprocess.check_output('gdal_translate -CO COMPRESS=deflate -of GTiff -a_srs EPSG:4326 -a_ullr '+str(xmin)+' '+str(ymin)+' '+str(xmax)+' '+str(ymax)+' -a_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" NETCDF:"'+filepath+'":Rw490 "'+filepathout+'Rw490_'+ID+'_10m.tif"', shell = "True") subprocess.check_output('gdalwarp -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" -t_srs "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"', shell = "True") f_in.close()
I attached an image to illustrate my problem. In the image, the red line represents the "good" position in lat/lon and it's superimposed to the Polymer result converted in Geotiff. I thought the problem was due to the product S2 but I tried on several products (different date and different location), and the problem is the same.
Did anyone already try to convert NETCDF to GeoTIFF or have this problem?