Geocoding S-1 images

[nukefix]

This thread is for discussing S-1 geocoding issues

[nukefix]

Here I've copied info about the required UTM-projection. What is the desired pixel spacing?

Yes. Landsat-8 uses Universal Transverse Mercator (UTM) map projection for everything except south of -63º latitude (ie Antarctica) where Polar Stereographic is used, both relative to World Geodetic System 84 datum (WGS84).

Petermann is UTM_ZONE = 20, less commonly 21 as per the metadata file ending in MTL.txt in each package download 
 
Jakobshavn is UTM_ZONE = 22, less commonly 23

Zachariae is UTM_ZONE = 27 for the area Wipneus and Espen have been looking at

Sptizbergen, Jan Mayen, and Iceland we don't look at too often.

[A-Team]

desired pixel spacing?

I would say each pixel should represent 15m x 15m on the ground, same as the nominal resolution of panchromatic band 8 of Landsat-8. We sometimes use other bands in the visible which are 30m but those are easily adjusted to 15m; indeed following Wipneus' protocol, we do that for 'pan-sharpening'.

Landsat provides corner and center lat,lon but no grid; typically we would crop severely to a region of interest and lose these known points. There is also a file available for download called "LandsatLook images with Geographic Reference (13.4 MB)" which sounds promising (am just downloading one now).

This folder contains various thermal infrared clouds, the 731 preview file, kml to place them on google earth, and a number of ".wld" files which are not images but "world" format (a plain text data file used by geographic information systems GIS to georeference raster map images, introduced by Esri) which does not describe the projection but is explained further at https://en.wikipedia.org/wiki/World_file and look like this opened in a text editor:

30.0000000000 pixel size in the x-direction in map units/pixel
0.0000000000 rotation about y-axis
0.0000000000 rotation about x-axis
-30.0000000000 pixel size in the y-direction in map units, almost always negative
417885.0000000000  x-coordinate of the center of the upper left pixel UTM northing
9099615.0000000000 y-coordinate of the center of the upper left pixel UTM easting

What I am getting at here is after getting Sentinel into UTM with 15m pixels, then we want to co-register it precisely with Landsat using geocoding, rather than manually align in Gimp etc., allowing for a script.

Just to cross-post a bit, here is some very promising information I received today from the help desk representative at Copernicus EO Support:

"To complete our previous message, the user can see web ref https://scientiaplusconscientia.wordpress.com/2015/08/06/working-with-sentinel-1-data-pre-processing-georeferencing-and-exporting-with-snap/
This even describes how to superimposes a Landsat and S1A image. The blog also describes the ESA tool SNAP (http://step.esa.int/main/toolboxes/snap/)."

Sentinel-1 imagery is provided in satellite geometry which is slightly different depending on the type of Level 1 product. For the detected GRD products, the imagery is provided in ground range, x, and azimuth, y, geometry whereas for complex SLC products, the imagery is provided in slant range, x, and azimuth, y, geometry. This is described on page https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/product-types-processing-levels/level-1.

As both GRD and SLC imagery are provided in satellite geometry, its orientation with respect to true north will vary with latitude and incidence angle. Thus it is not possible to accurately orientate the imagery with respect to north using a single rotation as this depends on both the x (range) and y (azimuth) position within the image.

Although Sentinel-1 imagery is not provided in a geocoded orientation (like your Landsat image), the product annotation gives a complete description of the transformation from pixel coordinates (x,y) to latitude & longitude via the geo-location data set record. This data set is described in the Sentinel-1 Product Specification document S1-RS-MDA-52-7441 available from https://sentinel.esa.int/documents/247904/349449/Sentinel-1+Product+Specification+3.0. In particular see Section 6.3.1.7 for a description of the geolocation grid. Below is an example of first point within a grid from a GRD product:

<geolocationGrid>
<geolocationGridPointList count="483">
<geolocationGridPoint>
<azimuthTime>2015-08-04T22:03:55.263139</azimuthTime>
<slantRangeTime>4.976224778237438e-03</slantRangeTime>
<line>0</line>
<pixel>0</pixel>
<latitude>7.768200950286953e+01</latitude>
<longitude>-9.014424336694947e+01</longitude>
<height>-2.077938625589013e-01</height>
<incidenceAngle>1.937838349149362e+01</incidenceAngle>
<elevationAngle>1.737100868906543e+01</elevationAngle>
</geolocationGridPoint>

Note that an approximate orientation of the image with respect to north can be found by extracting the platformHeading parameter within the product annotation (see Table 6-31 of the above product specification document). Note that this is the orientation of the satellite and not the image. An example of this parameter is given for an image at a high latitude:

<generalAnnotation>
<productInformation>
<pass>Ascending</pass>
<timelinessCategory>Fast-24h</timelinessCategory>
<platformHeading>-4.390998067621028e+01</platformHeading>

[nukefix]

The image I just posted in the Jakobshavn-thread should have the same geometry as the Landsat scenes. Of course, since orthorectification of the image was not possible (no DEM!) there are still altitude-dependent distortions in the SAR-image. This should not matter much at the calving-front though.

[nukefix]

The info from the helpdesk looks technically correct but it is somewhat beside the point - in order to do coordinate-transformation on SAR images one should use the right tools, for example SNAP (previously known as the S1TBX).

In SNAP one needs to:

1. Open downloaded .zip
2. Select 'Calibrate' under the Radar/Radiometric - menu item. Generating Sigma0 is enough.
3. Ellipsoid geocode the calibrated image by selecting 'Average-Height Range-Doppler' under Radar/Geometric/Ellipsoid Correction - menu item. Set the right projection and pixel-size in the dialog.
4. Add log-intensity to ellipsoid corrected calibrated image by selecting the Sigma0-band and selecting 'Linear to/from dB' under the Raster/Data Conversion - menu item.

This is rather simple and quick with a decent workstation. Of course, all of this can be scripted so that it can be run automatically from the command-line. If there was a high-resolution DEM available it would be equally simple to use it for orthorectifying the SAR-data.

Regarding superimposing Ellipsoid gecoded images with Landsat I'd expect that there will be some kind of a shift between the images. I'm guessing a single hand-picked GCP will do the trick.

[A-Team]

with a decent workstation

I read this as saying Copernicus should do the computation and provide an extra to link these images in geocoded Landsat compatible form as this is by far the most common thing for user to want and is far more efficiently computed centrally on a supercomputer, by purchasing cpu cycles and storage cheaply at Amazon AWS, or doing as a distributed project.

Alternatively -- and what has happened to Landsat -- is some enterprising but disgusted individual simply downloads all the imagery and sets up a separate site that provides the information in the form that users want. In the case of Landsat, users want to pick and choose band downloads, say band 8 only, and there only to the boundary coordinates of a rectangular mask. This is hugely more efficient from the end user perspective.

[nukefix]

In order to save processing effort in SNAP (the Sentinel toolboxes) it is possible to chain the required operations in a tool called Graph Builder and then execute them in one go. Here's a screenshot of a graph that calibrates S-1 images, geocodes them into UTM Zone 22 and converts the output from linear scale into decibels.

[A-Team]

Nice. Much more conducive for user engagement.

Reminds me of a tear-down someone asked me to do of Snowdon spyware.
http://electrospaces.blogspot.com/2014/03/olympia-how-canadas-csec-maps-phone-and.html

[nukefix]

Updated processing graph in SNAP. Run a similar chain on a bunch of images using batch-processing functionality and then coregister with create stack under radar - excellent results pretty much guaranteed.

[A-Team]

Run a similar chain on a bunch of images using batch-processing functionality and then coregister with create stack under radar ... excellent results pretty much guaranteed.

What is the time scale needed for batch processing say 20 images each for the main glaciers Pine Island, Thwaites, Totten, Zachariae, Petermann and Jakobshavn?

Does it let you crop the stack to a mask to save on file size or processing time? What are the options for saving a stack, I suppose these are 8-bit tif in the case of Sentinel 1A so perhaps can save as a basic gif animation.

It seems like one whole Sentinel scene is already fairly bulky at 10 m resolution, so a stack (and various derived products) is way too much for the forum. It sounds like we should set up cloud storage linked to by forum thumbnails rather than have each person repeat all the processing (which few would do). It is sort of an access issue ... how can people with limited desktop resources still participate?

I see they released SNAP 2.0 the other day. Did this fix the issue you were having with polar stereographic?

As your processing graphs become more complex, do they become portable in the sense of something others can copy/paste somewhere into SNAP, or do they have to replicate your graphs icon by icon?

The graphical user interface for the Mac is very well done, I've just opened a Landsat-8. This tool does a ton of things. ESA is going to "own" this whole business if they keep this up.

The Sentinel-1 Toolbox (S1TBX) consists of a collection of processing tools, data product readers and writers and a display and analysis application to support the large archive of data from ESA SAR missions including SENTINEL-1, ERS-1 & 2 and ENVISAT, as well as third party SAR data from ALOS PALSAR, TerraSAR-X, COSMO-SkyMed and RADARSAT-2.

The various processing tools could be run independently from the command-line and also integrated within the graphical user interface. The Toolbox includes tools for calibration, speckle filtering, coregistration, orthorectification, mosaicking, data conversion, polarimetry and interferometry.

The Sentinel-1 Toolbox is being developed for ESA by Array Systems Computing in partnership with DLR, Brockmann Consult and OceanDataLab.

http://step.esa.int/main/snap-2-0-out-now/

[nukefix]

On my kick-ass 6-core workstation processing a single scene takes about 2.5 minutes. I suppose I could subset with coordinates or shapefile but I have not tried that yet. It is possible to reduce the number of bits per pixel but I haven't checked the options out. Polar Stereo still does not work in the south. Graphs are straight XML so they can be passed on to others easily.

[A-Team]

processing a single scene takes about 2.5 minutes

And then it can go unattended in batch mode at this 24/hour rate while you are putzing with email or not even there? Decent utilization of a computer.

A shape file for certain glaciers might save on processing time but not so much for PI because the tributaries are so sprawled out. Masking the end product to the ROI might halve final file size given how png and jpg compress blocks.

In images that you posted like PIG_stack2_RGB1.jpg hair dresser, there was full use of 8-bit in each of the three channels (days). Running this through 2-bit, 3-bit, ... 7-bit, there is noticeable loss of image quality until 8-bit is hit. So no savings there. The binary has some use as a mask for enhancing  8-bit clarity.

I was just watching this very worthwhile E Rignot video on Pine Island from 2013, https://www.youtube.com/watch?v=2fMnauMEgxw&feature=youtu.be   He says two-thirds through that it is important and unexpected that the entire drainage basin up to the crest responds to drawdown at the coast -- yet they have been unable to demonstrate that yet for PI for lack of data.

So that would be very cool if we could show that using the earliest and latest Sentinels (which would maximize sensititivity).  However intermediate dates provide important internal consistency and quality controls on seasonality etc. There might be 15 months of data? Better if there were two full years.

Whatever, it should be possible to pick up very small amounts of motion if there is a decent level of distinctive features. This may require other orbital scenes than the one you have been processing which is slightly cut off even on the arches.

From #458, the overall scene image is PIG_stack2_RGB3.jpg is 760 KB, 2132x1269 pixels but at 1/5 the resolution of the two 10 m. So a whole scene would be 19=25*0.761 MB  at 10660x6345.

At the far right of the central trunk image, the displacement is 19 pixels or 190 m per 48 days (blue to red displacement) or ~4.0 m/day or 1.45 km/year. This is a decent velocity given the distance to the coast is about half the distance to the tributary crests.

How fast is the velocity falling off with distance from the calving front? Near the far left of the central trunk image (which is 17. km wide), the displacement is larger at 21.8 pixels. So the velocity is falling off by ~ 30 m per 15 km per 48 days or 15 m per km per year down from the 1450 km per year or 1% of the velocity is lost per km.

If a 1 pxl displacement over a year between two scenes could optimistically be detected, then the ultimate sensitivity is 76 m/yr.  For comparison, just below the crest at Summit, Greenland the velocity is 1 m/yr.

[Andreas Muenchow]

BIG 120 day BUMP:

Trying to follow this procedure, I run into trouble that the GeoTIFF does not display anywhere except within the software that generated it (ESA's SNAP Processor), even though gdalinfo executed on an X11 command line gives me

[wifi-roaming-128-4-223-79:~/Desktop] a2% gdalinfo S1B_output.tif
Driver: GTiff/GeoTIFF
Files: S1B_output.tif
Size is 14206, 13515
Coordinate System is:
PROJCS["WGS 84 / UTM zone 4N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-159],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32604"]]
Origin = (526.640284275286831,8193473.883620312437415)
Pixel Size = (40.000000000000000,-40.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_IMAGEDESCRIPTION=S1B_output
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (     526.640, 8193473.884) (174d41' 2.19"W, 73d14'33.10"N)
Lower Left  (     526.640, 7652873.884) (171d18' 4.29"W, 68d32'12.98"N)
Upper Right (  568766.640, 8193473.884) (156d47'18.09"W, 73d49'19.82"N)
Lower Right (  568766.640, 7652873.884) (157d16'55.59"W, 68d58'40.88"N)
Center      (  284646.640, 7923173.884) (165d 1'53.03"W, 71d18'50.62"N)
Band 1 Block=14206x13515 Type=UInt16, ColorInterp=Gray
[wifi-roaming-128-4-223-79:~/Desktop] a2%

which is "almost" identical to a LandSat GeoTIFF that does display correctly. The only difference between the working LandSat and the not-working Sentinel-1B is that "AREA-OR-POINT" is "AREA" for Sentinel-1B while it is "POINT" for LandSat ... and the dot-per-inch appears as 1 pixel/inch for Sentinel while it is 72 pixel/inch for LandSat.

Anyone know what I am doing wrong here? [Apologies if this is a newby question, perhaps, but I only got into this radar processing business 2 days ago.]

[Tealight]

I never used GeoTIFF but I know that SNAP also has an "Export" function under File. When I used it S1 images displayed correctly in GeoTIFF. Maybe it helps you to find the error and build a working graph.

[Andreas Muenchow]

Yes, I exported data as geotiff files. It turns out that my subsequent (gdal and fortran) codes have no problem with this format, but Apple's Preview has, if the files are too large. If I choose the subsetting to generate smaller geotiff files, Preview works again.

So, SNAP and Sentinel-1 are awesome products ... Preview not so much.