1

I am trying to convert Sentinel-3 SLSTR S6 radiance data to a Geotiff file. 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.

The code that I am trying before now is below.

Can somebody guide me in the process?

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)

2 Answers 2

0

When using Netcdf4 data for reprojection or assigning projection, using Pyresample library worked best for me. Especially kd_tree and geometry extension. You can also find a sample code snippet here. They use lat long info from another file just as you are doing and then reproject the file using Pyresample.

Here is a function that helped me in testing the library and understand how it works.I am using another tile of a similar area although a different dataset (with the right projection) to then reproject my current tile (without correct projection). Might not be directly useful, however can get you going with a few changes as per your requirement.

def proj_assign(proj_data,swath_lons,swath_lats,proj_info,width,height):
    area_id = 'test'
    description = 'test area'
    proj_id = 'test'
    #proj_info = assign_proj(tile)
    projection = get_projection(proj_info)
    lons = np.arange(proj_info[0][0],proj_info[0][0]+proj_info[0][1]*width,proj_info[0][1])
    lats = np.arange(proj_info[0][3],proj_info[0][3]+proj_info[0][5]*height,proj_info[0][5])
    area_extent = (lons[0], lats[-1], lons[-1], lats[0])
    area_def = geometry.AreaDefinition(area_id, description, proj_id, projection, width, height, area_extent)
    swath_def = geometry.SwathDefinition(lons=swath_lons, lats=swath_lats)
    result = kd_tree.resample_nearest(swath_def, proj_data, area_def)
    return result, proj_info
-1

You could use high level libraries such as EOReader to geocode directly your data:

from eoreader.reader import Reader
from eoreader.bands import SWIR_2 # which is S6 band

# Open the Sentinel-3 product
s3_path = "S3B_SL_1_RBT____20191115T233722_20191115T234022_20191117T031722_0179_032_144_3420_LN2_O_NT_003.SEN3"

prod = Reader().open(s3_path)

# Load S6 band
s6 = prod.load(SWIR_2)[SWIR_2]
s6.plot()

enter image description here

More information on Sentinel-3 data are available in this notebook.

Disclaimer: I am the maintainer of EOReader

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.