1

Let's assume a netCDF4 (.nc) file downloaded from https://s5phub.copernicus.eu/dhus/odata/v1/Products('e6e91b26-ca43-43d4-9c08-c9c14dd6737e')/$value . This is sulphur dioxide measurement dataset over Vietnam, China and Filippines, given in the linked image.

Sulphur dioxide data over Vietnam, China and Filippines

I wish to transform this into a GeoTIFF file for further processing so that I can process the values in (lon, lat, value) manner. I start by the following command line gdal_translate -ot float32 -unscale -CO COMPRESS=deflate HDF5:"S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719.nc"://PRODUCT/sulfurdioxide_total_vertical_column so2.tif.

When I open the resulting GeoTIFF file, it's almost pitch black with some little grey dots. Using Panoply to plot a diagram produces the following

enter image description here

I also see in the metadata there reads

float sulfurdioxide_total_vertical_column(time=1, scanline=278, ground_pixel=450);
:_FillValue = 9.96921E36f; // float
:units = "mol m-2";
:standard_name = "atmosphere_mole_content_of_sulfur_dioxide";
:long_name = "total vertical column of sulfur dioxide for the polluted scenario derived from the total slant column";
:coordinates = "/PRODUCT/longitude /PRODUCT/latitude";
:multiplication_factor_to_convert_to_DU = 2241.15f; // float
:multiplication_factor_to_convert_to_molecules_percm2 = 6.02214E19f; // float
:_ChunkSizes = 1U, 278U, 450U; // uint

Question:

Do I make this right so that the data values are correct if I intend to use them later as (lon, lat, value) (with WGS84 projection, read from metadata) in various calculations?

E.g. that I don't accidentally corrupt the coordinates or the data values (by accidentally scaling, say). Put otherway, I'm a n00b and asking if this method of using the tooling is OK before banging my head further. I'm still in process of learning the other parts of data processing and checking the correctness of those steps would be easier if I were sure about this step.

<edit 1: Looking at the tiff metadata, it doesn't appear the data is georeferenced. However gdalwarp so2.tif so2georef.tif -t_srs "+proj=longlat +ellps=WGS84" fails with

ERROR 1: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for so2georef.tif. There is no affine transformation and no GCPs. Specify transformation option SRC_METH OD=NO_GEOTRANSFORM to bypass this check.

It think this is further confirmed by gdal_translate -of XYZ so2.tif so2.xyz that shows values beginning with

0.5 0.5 9.96920996838686905e+36
1.5 0.5 9.96920996838686905e+36
2.5 0.5 -2.67918767349328846e-06
3.5 0.5 0.000152658714796416461
4.5 0.5 -0.000818965549115091562
5.5 0.5 0.000193428524653427303
6.5 0.5 0.00010718178964452818
7.5 0.5 -0.00120470335241407156
8.5 0.5 0.000939452263992279768

So either something is wrong with the initial command or the original .nc file data needs to be augmented. I'm currently too rookie to understand all the metadata but I do see the portal has shapefiles etc.

<edit 2: It appears the glitch was about environment variables and doing a bit of fiddling with that, a new command line I tried is gdal_translate -ot float32 -unscale -CO COMPRESS=deflate -of GTiff -a_srs EPSG:4326 HDF5:"S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719.nc"://PRODUCT/sulfurdioxide_total_vertical_column so2.tif

So I added -of GTiff -a_srs EPSG:4326 and it appears the values I get to the xyz file are still the same. It appears the metadata now tells the coordinates are georeferenced, but they aren't really transformed as should, just the metadata added without real transformation. In the metadata I have spotted

geolocation_grid_from_band=3
geospatial_lat_max=25.532017
geospatial_lat_min=2.1418159
geospatial_lon_max=99.929924
geospatial_lon_min=128.72141

This is probably data that can be explicitly used. If so, might be nice to know if gdal_translate could automatically use it, perhaps with an explicit switch.

<edit 3: There's S-P instruction manual at http://www.tropomi.eu/sites/default/files/files/S5P-KNMI-L2-0009-SD-S5P_level_2_Input_Output_Data_Definition-9.0.0-20170614.pdf, from there the following explanation. Then judging from https://www.gdal.org/gdal_translate.html maybe -a_srs EPSG:4326 -a_ullr 2.1418159 128.72141 25.532017 99.929924 could be added (or in some other order) to the original command line.

S-5P coordinate referencing explanation

  • 1
    Welcome to GIS SE. As a new user, please take the Tour, which emphasizes the importance of asking one question per Question. Please Edit this question to focus on one answerable question. – Vince Oct 31 '18 at 2:45
  • Check! I removed the more nice-to-know-questions and concentrated only to one question with a reason to ask it. – Veksi Oct 31 '18 at 5:29
1

It appears the approximately command line to do this is gdal_translate -ot float32 -unscale -CO COMPRESS=deflate -of GTiff -a_srs EPSG:4326 -a_ullr 25.532017 99.929924 2.1418159 128.72141 HDF5:"S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719.nc"://PRODUCT/sulfurdioxide_total_vertical_column so2.tif

where the a_ullr parameters are derived from the .nc file metadata, reproduced here:

geolocation_grid_from_band=3
geospatial_lat_max=25.532017
geospatial_lat_min=2.1418159
geospatial_lon_max=99.929924
geospatial_lon_min=128.72141

Applying these produces a geotiff file with the following georeferecing information

Corner Coordinates:
Upper Left  (      25.532,      99.930) ( 25d31'55.26"E, 99d55'47.73"N)
Lower Left  (      25.532,     128.721) ( 25d31'55.26"E,128d43'17.08"N)
Upper Right (       2.142,      99.930) (  2d 8'30.54"E, 99d55'47.73"N)
Lower Right (       2.142,     128.721) (  2d 8'30.54"E,128d43'17.08"N)
Center      (      13.837,     114.326) ( 13d50'12.90"E,114d19'32.40"N)

Where the Upper Left and Lower Right are approximately correct but otherwise Lower Left and Upper Right have switched places. Also comparing the ESA satellite image (linked) and then inserting the corner coordinates to Google Maps yield slightly different places on the map. I tried to project with EPSG:3857 already and I tried placing the corder coordinates in different order, no dice. That's a different set of questions though. I leave this here in case someone finds it useful, a follow up at GeoTIFF (ESA Copernicus data) georeferencing coordinates are mirrored and projection slightly off when checking on Google Maps .

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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