4

I worked on downloading a Sentinel 5P data in my question at Converting a netCDF4 file to (georeferenced) GeoTIFF, problems with georeferencing and also transforming the netCDF4 (.nc) file into a geoTIFF image/raster data set.

As noted in the answer to my question, it appears the georeferencing method I use doesn't work correctly.

  1. Two of the corner points are flipped and for a reason or another I don't seem to get them correctly.
  2. I use Google Maps to check I have transformed the geoTIFF file correctly. In addition. The corner points are also slightly off of what Copernicus map shows. I assume this is a difference in the way they project the maps, but as the ultimate goal is to get the tiff on a common web map, I think I need to get this part right too.

Question: Can someone shed light on what's wrong with my projection parameters? The parameters are given in the following, first a bit of setup to make it easier to undersand what I was thinking.

I replicate the previous question here as appropriate for this one. In case it helps, the Copernicus portal is at https://s5phub.copernicus.eu/dhus/#/home. The guest username and password are given in a prompt when one presses login. It has a lot more tooling to fiddle with ideas, shapefiles amongst other things.

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, Taiwan and Philippines, given in the linked image. Regarding the portal, searching with S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719 yields the same data/file as downloading directly.

Sulphur dioxide data over Vietnam, China and Filippines

From the downloaded .nc file I can see in the metadata

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

so I proceed the create a geoTIFF raster data file as follows: 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

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 (as an aside, it'd be nice to know if this can be done without dumping metadata about the corner points first and applying the information specifically).

Here are the Google Maps links for the corner points for convenience:

<edit: With the proposed parameters -a_ullr 99.930 25.532017 128.72141 2.142 gdalinfo so2.tif gives

Corner Coordinates:
Upper Left  (  99.9300000,  25.5320170)
Lower Left  (  99.9300000,   2.1420000)
Upper Right ( 128.7214100,  25.5320170)
Lower Right ( 128.7214100,   2.1420000)
Center      ( 114.3257050,  13.8370085)

which doesn't seem be right. Looking at the coordinates, it gives a feeling the corners are on straight lines, like a square on a flat map. I'll be spelunking something with "spherical Mercator" soon, I suppose. Basically waving hands to every direction and seeing if something sticks. :P

<edit 5: Maybe it's something about snapping values to pixel colors, i.e. a drawing artifact. Here's another way of seeing the image after removing sentinel values:

SO2 image, without the black band

<edit 4:

On QGIS looking at metadata about sentinel values and setting their opacity to 100% crops of the edges and reveals a black image that looks like the original. Then adjusting the histogram so that the maximum value is the highest measurement point gives the grey scale image in the picture. The black band on the left side seems to be all minimum values. I don't know the cause of this. It could be a projection artefact but it could also be invalid QA values that should be cropped off. Here's the image:

sentinel value opacity set to 100% and histogram adjusted

<edit 3:

I proceed with the instructions hinted by @AndreJ and run the following commands

gdal_translate -of VRT HDF5:"S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719.nc"://PRODUCT/longitude lon.vrt
gdal_translate -of VRT HDF5:"S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719.nc"://PRODUCT/latitude lat.vrt
gdal_translate -of VRT HDF5:"S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719.nc"://PRODUCT/sulfurdioxide_total_vertical_column so2.vrt

And modify the resulting so2.vrt to read as follows

<VRTDataset rasterXSize="450" rasterYSize="278">
   <Metadata>
   <!-- Snipped off for brevity. -->
   </Metadata>
   <!-- Added this geolocation domain information. -->
   <Metadata domain="GEOLOCATION">
       <mdi key="X_DATASET">lon.vrt</mdi>  
       <mdi key="X_BAND">1</mdi>  
       <mdi key="Y_DATASET">lat.vrt</mdi>  
       <mdi key="Y_BAND">1</mdi>  
       <mdi key="PIXEL_OFFSET">0</mdi>  
       <mdi key="LINE_OFFSET">0</mdi>  
       <mdi key="PIXEL_STEP">1</mdi>  
       <mdi key="LINE_STEP">1</mdi>
   </Metadata>    
   <VRTRasterBand dataType="Float32" band="1">
   <Metadata>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column_coordinates">/PRODUCT/longitude /PRODUCT/latitude</MDI>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column_long_name">total vertical column of sulfur dioxide for the polluted scenario derived from the total slant column</MDI>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column_multiplication_factor_to_convert_to_DU">2241.1499 </MDI>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column_multiplication_factor_to_convert_to_molecules_percm2">6.02214e+19 </MDI>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column_standard_name">atmosphere_mole_content_of_sulfur_dioxide</MDI>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column_units">mol m-2</MDI>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column__FillValue">9.96921e+36 </MDI>
       <MDI key="PRODUCT_sulfurdioxide_total_vertical_column__Netcdf4Dimid">2 </MDI>
     </Metadata>    
     <SimpleSource> 
         <SourceFilename relativeToVRT="1">HDF5:S5P_NRTI_L2__SO2____20181030T054006_20181030T054506_05418_01_010102_20181030T062719.nc://PRODUCT/sulfurdioxide_total_vertical_column</SourceFilename>
         <SourceBand>1</SourceBand>
         <SourceProperties RasterXSize="450" RasterYSize="278" DataType="Float32" BlockXSize="450" BlockYSize="278" />
         <SrcRect xOff="0" yOff="0" xSize="450" ySize="278" />
         <DstRect xOff="0" yOff="0" xSize="450" ySize="278" />
      </SimpleSource>
 </VRTRasterBand>
 </VRTDataset>

Then running gdalwarp -geoloc -t_srs EPSG:4326 -overwrite so2.vrt so2.tif and finally gdalinfo so2.tif yields

Corner Coordinates:
Upper Left  ( 100.1267700,  25.0621347) (100d 7'36.37"E, 25d 3'43.69"N)
Lower Left  ( 100.1267700,   2.5234540) (100d 7'36.37"E,  2d31'24.43"N)
Upper Right ( 128.5601826,  25.0621347) (128d33'36.66"E, 25d 3'43.69"N)
Lower Right ( 128.5601826,   2.5234540) (128d33'36.66"E,  2d31'24.43"N)
Center      ( 114.3434763,  13.7927944) (114d20'36.51"E, 13d47'34.06"N)
Band 1 Block=493x4 Type=Float32, ColorInterp=Gray

which isn't correct. Interestingly the coordinates are "inverted" in that changing 100.1267700, 25.0621347 -> 25.0621347, 100.1267700 is approximately correct, albeit slightly offset when compared to the Copernicus image (looking the coordinates via Google Maps now). It doesn't appear the reason is solely about longitude and latitude information changed but also the projection is wrong.

Opened in QGIS and projected on Google Maps gives SO2 geotiff file which looks like being in the right place. The data is unscaled so nearly black (I assume), the white strand is an interesting looking thing (looks like being the edge of the original satellite picture, i.e. no data, but why black on the other side then?).

<edit 2: There's interesting discussion at https://github.com/stcorp/harp/issues/189 about possible issues. Even if the discussed issue isn't directly this, one should note it's about Sentinel 5P, georeferencing and invalid values in certain situations. For that I wonder if that's the case here and how to find those values and filter out. As a tangential educative issue, pay attention to vocabulary about level 1/2/3 products too, see more at https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-5p/products-algorithms.

<edit 1: In the metadata I see (I don't know if this matters):

PRODUCT_corner_CLASS=DIMENSION_SCALE
PRODUCT_corner_comment=This coordinate variable defines the indices for the pixel corners; index starts a 0 (counter-clockwise, starting from south-western corner of the pixel in ascending part of the orbit).
PRODUCT_corner_long_name=pixel corner index
PRODUCT_corner_NAME=corner
PRODUCT_corner_REFERENCE_LIST=
PRODUCT_corner_units=1
PRODUCT_corner__FillValue=-2147483647
PRODUCT_corner__Netcdf4Dimid=3
PRODUCT_delta_time_long_name=offset from reference start time of measurement
PRODUCT_delta_time_units=milliseconds
PRODUCT_delta_time__FillValue=-2147483647
PRODUCT_delta_time__Netcdf4Dimid=2
PRODUCT_ground_pixel_axis=X
PRODUCT_ground_pixel_CLASS=DIMENSION_SCALE
PRODUCT_ground_pixel_comment=This coordinate variable defines the indices across track, from west to east; index starts at 0
PRODUCT_ground_pixel_long_name=across-track dimension index
PRODUCT_ground_pixel_NAME=ground_pixel
PRODUCT_ground_pixel_REFERENCE_LIST=
PRODUCT_ground_pixel_units=1
PRODUCT_ground_pixel__FillValue=-2147483647
PRODUCT_ground_pixel__Netcdf4Dimid=1
PRODUCT_latitude_bounds=/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/latitude_bounds
PRODUCT_latitude_long_name=pixel center latitude
PRODUCT_latitude_standard_name=latitude
PRODUCT_latitude_units=degrees_north
PRODUCT_latitude_valid_max=90
PRODUCT_latitude_valid_min=-90
PRODUCT_latitude__FillValue=9.96921e+36
PRODUCT_latitude__Netcdf4Dimid=2
PRODUCT_layer_CLASS=DIMENSION_SCALE
PRODUCT_layer_long_name=layer dimension index
PRODUCT_layer_NAME=layer
PRODUCT_layer_REFERENCE_LIST=
PRODUCT_layer_units=1
PRODUCT_layer__Netcdf4Dimid=4
PRODUCT_longitude_bounds=/PRODUCT/SUPPORT_DATA/GEOLOCATIONS/longitude_bounds
PRODUCT_longitude_long_name=pixel center longitude
PRODUCT_longitude_standard_name=longitude
PRODUCT_longitude_units=degrees_east
PRODUCT_longitude_valid_max=180
PRODUCT_longitude_valid_min=-180
PRODUCT_longitude__FillValue=9.96921e+36
PRODUCT_longitude__Netcdf4Dimid=2
  • 2
    GDAL handles coordinates in longitude-latitude and easting-northing order, correct that. But it is odd that in the metadata lon_min is bigger than lon_max. – user30184 Oct 31 '18 at 13:21
  • 1
    Think "upper-left = least east, most north". So -a_ullr 99.930 25.532017 128.72141 2.142. But it can be that your image can't be georeferenced properly with just two corner coordinates but you should see it soon. – user30184 Oct 31 '18 at 13:45
  • 1
    It looks like your data has longitude and latitude bands to use. Since the image is bended on your map, -a_ullr is not the right tool. You can use it only if the image is expected to be parallel to longitudes and latitudes. – AndreJ Nov 2 '18 at 10:39
  • 1
  • 1
    You might as well look into my answers at gis.stackexchange.com/questions/154339/… and gis.stackexchange.com/questions/81361/… – AndreJ Nov 3 '18 at 15:37
1

I think you are already there. I got the same extent as you, and this image:

enter image description here

Note that the extent of the warped tif is now a rectangle in EPSG:4326, with corners different to the bended satellite view.

The blue dashed line I got from the METADATA_EOP_METADATA_om:featureOfInterest_eop:multiExtentOf_gml:surfaceMembers_gml:exterior_gml:posList variable (which uses lat-lon order).

Keep in mind that cell center coordinates and raster extent differ by half the cell size.

  • So that blue dashed line is in WGS84 system to begin with but I just didn't realize it (since the metadata is in .nc)? :) Then if I want to show it on Google Maps, or whatever, then I do the translation and show the square image with, say, alpha channel for edges where there really isn't data. That is an excellent note about the square and the different order. I could actually accept this an answer, but if I may ask a final question, shouldn't I produce a geotiff image with the coordinates inverted so that Upper Left etc. are correct when checked on maps? – Veksi Nov 4 '18 at 15:39
  • I added a further comment and a image. I don't understand why the white strand in the image. It looks like it goes on the edge of actual measurement data, but why is its left side (the western side) then again black? It feels like there's something wrong in the transformation still, but I don't know. It feels strange it's all black if there isn't measurement data. I think I need to check the fill values or something to clarify this... – Veksi Nov 4 '18 at 16:42
  • 2
    I used -srcnodata 9.96921e+36 -dstnodata 9999 in my command line. – AndreJ Nov 4 '18 at 19:52

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.