CryoTEMPO SI.9/PO.9 - Implement a more accurate land-ocean mask

[shendric]

The current land/ocean mask is based on GSHGG data, which is incorrect. The task includes the following actions:

• Identify a suitable land/ocean mask
• Implement replacement of surface type flag in l1 pre-processor

[shendric]

Which landmask and type (grid/shapefile) should we use?

[shendric]

This issue also requires the change of the Level-1 pre-processor:

At the moment there is no injection point of auxiliary information, before the extraction of polar ocean segments. This method already relies on the existing land/ocean flag in the l1 source data and thus an injection point for data must be established directly after the source file is parsed.

This will also be beneficial if data should be added to the l1p data object that exists for the same period as the source data.

[shendric]

Three available land and or ocean shapefile products have been downloaded:

1. Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) v2.3.7
2. Natural Earth (NE) 1:10 million coastline
3. Open Street Map (OSM) Land Polygons

The resolution of the GSHHG (at full resolution) and OSM exceed the NE coastline. However a critical test is the coast line of greenland and all coastlines were evaluated against Global 30 Arc-Second Elevation (GTOPO30) digital elevation model:

The result was that the Natural Earth land mask hat the best agreement with the topographic data. But since Natural Earth lacks detail (small islands) in other parts of the Arctic Ocean, the goal for the moment is to create a land/ocean mask based on OSM, except for Greenland, which is based on NE.

[shendric]

As a first step, Greenland is manually removed from the OSM shapefile with QGIS. The remaining shapes are then rasterized using gdal on a 0.005 x 0.005 degree global grid (approximately 100m at 80N):

gdal_rasterize -l land_polygons -burn 255 -ot Byte -of gtiff -tr 0.005 0.005 land_polygons.shp land_polygons_raster_osm_wgs84.tiff

The chosen resolution preserved small islands (example here is Hinlopenstretet in Svalbard):

[shendric]

The total filesize of the global geotiff (date type Byte) is 2.4GB, which is excessive to hold in memory during the Level-1 pre-processor. The global data is therefore reprojecedt to a polar projection (WGS 84 / NSIDC EASE-Grid 2.0 North, EPSG:6931) again using gdal:

gdalwarp -s_srs WGS84 -t_srs EPSG:6931 -tr 250 250 -te -5000000 -5000000 5000000 5000000 -multi -r near -ot Byte -overwrite land_polygons_raster_osm_wgs84.tiff land_polygons_raster_osm_epsg6931.tiff

The output is a 1.2GB geoTIFF file with 40000P x 40000L pixels, which includes the full CryoSat-2 subset (north of 45N) used for CryoTEMPO:

[shendric]

Workflow to merge the Natural Earth Greenland polygon into the Open Street Map shapefile:

1. Select and remove the Greenland layer from the OSM shapefile using QGIS
2. Select Greenland layer from NE shapefile and remove all other shapes using QGIS (Note: In the NE land shapefile all land is one layer and Greenland itself cannot be easily selected. That is possible for the coastlines though and after selecting only Greenland the line object is converted into a polygon)
3. Import the NE Greenland polygon to the edited OSM shapefile
4. Merge the layers and save as new shape file.
5. Repeats rasterize and reprojection steps from above

The result:

[shendric]

Gdal also offers a way to compute the distance to the coast with gdal_proximity.py. This might be useful addition for filtering in later evolutions.

gdal_proximity.py land_polygons_raster_osm_epsg6931.tiff distance_to_coast_raster_osm_epsg6931.tiff -ot Float32 -values 255 -distunits GEO

The result for full mask extent:

and a subset at full resolution in the Canadian Archipelago:

[shendric]

Updated the Level-1 pre-processor workflow in shendric@d488118. It is now possible to the specify a stage for a L1 pre-processor item to run:

Possible pre-processor stages:

1. post_source_file
2. post_ocean_segment_extraction
3. post_merge

The change was necessary, because the updated land/ocean flag needs to be swapped in before the ocean segment extraction. And before were called only after the ocean segment extraction.

[shendric]

Running first tests in the Level-1 Pre-Processor:

• Memory footprint is a few < 5 GB.
• Extraction of both variables per source file is approximately 0.005 seconds
• Minor differences to ESA land/ocean flag observed, with potential improvements in one test case:

Next step is larger test.

[shendric]

Feature is on operation. This issue will be closed for now, but may be re-opened for land/ocean mask udpates