I have experience working with Landsat, Modis and Sentinel-2 data; Right now I need to use Sentinel-3 data also to fuse with Sentinel-2 data, however, I am finding converting the Sentinel-3 data which is in NCDF4 format to Geotiff format correctly very difficult. I am trying to use snappy API from Python. Since the lat-long information is provided in a separate file, my main problem is how to overlay geo-coordinates and radiance values from S6 and write it in a tiff file with the correct projection and resolution. Any help would be highly appreciated. Can someone please point me to some relevant examples
I would suggest that you don’t try to read the NetCDF files separately but read the xfdumanifest.xml file with snappy.
This way you get a product The handling of the data is probably easier then.
Hi , Many thank for replying. I have been reading the product .SEN3 and extracting as bands the individual files .The code I am trying to extract the Tiff file is below. So should I read the .xml file in the readProduct part? Can you please help me understand where am I going wrong here. My final aim is to converting the S6 radiance data to a Tiff file
p=ProductIO.readProduct(‘S3B_SL_1_RBT____20181201T035905_20181201T040205_20181214T172444_0180_019_161_2520_LN2_O_NT_003.SEN3’)
rad13 = p.getBand(‘S6_radiance_ao’)
w = rad13.getRasterWidth()
h = rad13.getRasterHeight()
rad13_data = np.zeros(w * h, np.float32) # numpy array of longitude
rad13.readPixels(0, 0, w, h, rad13_data)
lat = p.getBand(‘x_ao’)
long = p.getBand(‘y_ao’)
w1 = lat.getRasterWidth()
h1=long.getRasterHeight()
lat_data = np.zeros(w1 * h1, np.float32) # numpy array of latitude
lat.readPixels(0, 0, w, h, lat_data)
long_data = np.zeros(w1 * h1, np.float32)
long.readPixels(0, 0, w, h, long_data) # numpy array of longitude
After I got the separate numpy arrays of lat-long and radiance values I have been trying to convert it into a geotiff file
xmin,ymin,xmax,ymax = [long_data.min(),lat_data.min(),long_data.max(),lat_data.max()
nrows,ncols = np.shape(rad13_data)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)
output_raster = gdal.GetDriverByName(‘GTiff’).Create(‘myraster.tif’,ncols, nrows, 1 ,gdal.GDT_Float32)
output_raster.SetGeoTransform(geotransform)
####I am getting a 0 value here so not able to proceed
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
###getting a 0 value here
output_raster.SetProjection( srs.ExportToWkt())
output_raster.GetRasterBand(1).WriteArray(rad_array)