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[1] #netcdf file path
filepathout = sys.argv[2] #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?
Thank you!
Sonia