Convert Sentinel-2 results NETCDF to GEOTIFF

Post Reply
samaaouch
Posts: 4
Joined: Fri Jun 05, 2020 1:32 pm
company / institution: ARGANS Ltd
Location: Brest

Convert Sentinel-2 results NETCDF to GEOTIFF

Post by samaaouch »

Hi,

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()
My program generates a GeoTiff georeferenced BUT the output has an offset of few degrees.
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
You do not have the required permissions to view the files attached to this post.
User avatar
fsteinmetz
Site Admin
Posts: 140
Joined: Fri Sep 07, 2018 1:34 pm
company / institution: Hygeos
Location: Lille, France
Contact:

Re: Convert Sentinel-2 results NETCDF to GEOTIFF

Post by fsteinmetz »

Dear Sonia,

It seems indeed that there is a mismatch between the georeferencing you have defined in your geotiff file, and the "good" coastline, which I understand comes from your GIS software.
But did you actually use anything related to georeferencing from the polymer output ? The variables lat and lon are not used, so it is not clear to me how you defined the georeferencing of your geotiff file, and how it is related to polymer's output. I suspect there is an issue in the commands gdal_translate or gdalwrap but I am not a proficient user of gdal so I can not give you advice on these commands.

Cheers,
François
samaaouch
Posts: 4
Joined: Fri Jun 05, 2020 1:32 pm
company / institution: ARGANS Ltd
Location: Brest

Re: Convert Sentinel-2 results NETCDF to GEOTIFF

Post by samaaouch »

Dear François,


Yes, you're right I forgot to share the line which used Polymer lon/lat, so this the commands that I used to georeference my tiff with the missing line:

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]
filepathout = sys.argv[2]
f_in = Dataset(filepath, 'r')

for v in f_in.variables :
  print(v, np.shape(f_in.variables[v][:]))

lon = f_in.variables['longitude'][:]  
lat = f_in.variables['latitude'][:]
Rw490 = f_in.variables['Rw490'][:]

xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]

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()
Thank you for your help, I will continue to search!

Regards

Sonia
Post Reply