Port of Antwerp from space perspective

The Hackathon is just around the corner and so is the release of Earth Observation Analysis Service (EOaaS) coming with HANA 2 platform.

This is a great opportunity to give an overview of EOaaS from a differnt perspective.

So before you start, make sure you have a YaaS account ready and got your authorization token. Further help on this topic can be found in this blog post. The examples in this post contain the string `ACCESS_TOKEN`, which needs to be replaced by your current authorization token.

Notes:

• For the time being EOaaS is in beta.
• The used feature cloudPercentage is not available in beta and thus not yet documented in API doc.
• The format of CoverageId content is subject to change and will vary in future. However its functioinality won’t change.
• We use curl in a Linux shell. But feel free to use other tools.
• Currently only XML output is supported. In the near future, EOaaS will also support JSON. We won’t look too deep into the the XML file structure. However if you plan to write a web application, it might be worth waiting for the convenience coming with JSON.

Where am I and where do I get coordinates?

The ESA Copernicus program provides images all over the world. To retrieve only images from a certain region, EOaaS allows you to define a bounding box along latitude and longitude.

But before an appropriate bouding box can be created, we need to know where we are or from where we want images and therefore we need coordinates such as latitude and longitude.

In this blog post, we take a deeper look on the port of Antwerp. This is great, because we can find Antwerp on Wikipedia.

In this approach we get coordinates of Antwerp in degrees and minutes which is  51°13′N 04°24′E.

Alternatively a more accurate approach would be to use a web map interface to get the coordinates of a click location. An OpenLayers Tutorial on overlays does a
very good job. We zoom to Antwerp port and click in NW and SE corners, which reveal coordinates suites for an accurate bounding box.

We see the coordinates:

north west: 51°22′07"N 04°09′16"E
south east: 51°12′35"N 04°31′60"E

Note:
In both cases the coordinates are displayed in degree-minutes-second format. However we need decimal value. In our usecase it is sufficient to only to divide the minutes value by 60.

Simple conversion:
    51°22′ -> 51.3666  (51 remains, 22 is divided by 60)
Accurate conversion:
north west: 51.3686 N 4.1544 E
south east: 51.2097 N 4.5253 E

In case you have an address, such as “Waagnatie, hall 2 – Antwerp, Belgium” we can use a geocoder.

There are several commercial products in market. For now we try the open source alternative Nominatim and enter the address.

Unfortunately it does not find a result. But when we reduce the address string to “Waagnatie, Antwerp” we get a result. Now we open the map border window in top right corner and find the coordinates 51.23354, 4.40208.

Query images with boudning box

In the API documentation we find the parameter dimensionTrim which allows us to define the bounding box.

Using the upper retrieved value, we can attach the following paramters to our query:

dimensionTrim=lat(51.2,51.4)&dimensionTrim=long(4.1,4.6)

In our example we want standard RGB images, so we need to specify DS_RGB as
eoid`. This results in the following query URL.

https://api.yaas.io/sap/earth-analysis/v1/wcs/describeEoCoverageSet?eoid=DS_RGB&dimensionTrim=lat(51.2,51.4)&dimensionTrim=long(4.1,4.6)

So lets try it out! Get an authorization token (here called ACCESS_TOKEN) and run curl.

curl -H "Authorization: Bearer ACCESS_TOKEN" -o /tmp/Antwerp_query_result.xml -X GET "https://api.yaas.io/sap/earth-analysis/v1/wcs/describeEoCoverageSet?eoid=DS_RGB&dimensionTrim=lat(51.2,51.4)&dimensionTrim=long(4.1,4.6)"

We receive a long XML file with many wcs:CoverageDescription tag blocks.

In the moment of writing the first block looks as follows.

    <wcs:CoverageDescription gml:id="RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES">

        <gml:boundedBy>
            <gml:Envelope axisLabels="lat long"
                srsName="http://www.opengis.net/def/crs/EPSG/0/4326"
                uomLabels="deg deg"
                srsDimension="2">
                <gml:lowerCorner>50.4637984013 2.99971821167</gml:lowerCorner>
                <gml:upperCorner>51.4405411653 4.57954402146</gml:upperCorner>
            </gml:Envelope>
        </gml:boundedBy>

        <wcs:CoverageId>RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES</wcs:CoverageId>

        <gmlcov:metadata>
            <gmlcov:Extension>
                <wcseo:EOMetadata>
                    <eop:EarthObservation gml:id="eo_RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES">

                        <om:phenomenonTime>
                            <gml:TimePeriod gml:id="tp_RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES">
                                <gml:beginPosition>2016-11-24T10:43:52.000000</gml:beginPosition>
                                <gml:endPosition>2016-11-24T10:43:52.000000</gml:endPosition>
                            </gml:TimePeriod>
                        </om:phenomenonTime>

                        <om:resultTime />

                        <om:procedure>
                            <eop:EarthObservationEquipment gml:id="equ_RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES">
                                <eop:platform>
                                    <eop:Platform>
                                        <eop:shortName>SENTINEL2</eop:shortName>
                                    </eop:Platform>
                                </eop:platform>
                                <eop:instrument>
                                    <eop:Instrument>
                                        <eop:shortName>MSI</eop:shortName>
                                    </eop:Instrument>
                                </eop:instrument>
                                <eop:sensor>
                                    <eop:Sensor>
                                        <eop:sensorType>OPTICAL</eop:sensorType>
                                    </eop:Sensor>
                                </eop:sensor>
                            </eop:EarthObservationEquipment>
                        </om:procedure>

                        <om:observedProperty />

                        <om:featureOfInterest>
                            <eop:Footprint gml:id="footprint_RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES">
                                <eop:multiExtentOf>
                                    <gml:MultiSurface gml:id="ms_RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES" srsName="EPSG:4326">
                                        <gml:surfaceMember>
                                            <gml:Polygon gml:id="poly_RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES">
                                                <gml:exterior>
                                                    <gml:LinearRing>
                                                        <gml:posList>2.999712 51.451182 2.999718 50.463798 4.546436 50.453523 4.579544 51.440541 2.999712 51.451182</gml:posList>
                                                    </gml:LinearRing>
                                                </gml:exterior>
                                            </gml:Polygon>
                                        </gml:surfaceMember>
                                    </gml:MultiSurface>
                                </eop:multiExtentOf>
                            </eop:Footprint>
                        </om:featureOfInterest>

                        <om:result>
                            <cloudPercentage uom="%">100.0</cloudPercentage>
                        </om:result>

                        <eop:metaDataProperty>
                            <eop:EarthObservationMetaData>
                                <eop:identifier>RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES</eop:identifier>
                                <eop:acquisitionType>OTHER</eop:acquisitionType>
                                <eop:status>POTENTIAL</eop:status>
                            </eop:EarthObservationMetaData>
                        </eop:metaDataProperty>

                    </eop:EarthObservation>
                </wcseo:EOMetadata>
            </gmlcov:Extension>
        </gmlcov:metadata>

        <gml:domainSet>
            <gml:RectifiedGrid gml:id="RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES_grid" dimension="2">
                <gml:limits>
                    <gml:GridEnvelope>
                        <gml:low>0 0</gml:low>
                        <gml:high>10980 10980</gml:high>
                    </gml:GridEnvelope>
                </gml:limits>
                <gml:axisLabels>x y</gml:axisLabels>
                <gml:origin>
                    <gml:Point srsName="http://www.opengis.net/def/crs/EPSG/0/32631" gml:id="RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES_grid_origin">
                        <gml:pos>499980 5700000</gml:pos>
                    </gml:Point>
                </gml:origin>
                <gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/32631">10 0</gml:offsetVector>
                <gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/32631">0 -10</gml:offsetVector>
            </gml:RectifiedGrid>
        </gml:domainSet>

        <gmlcov:rangeType />

        <wcs:ServiceParameters>
            <wcs:CoverageSubtype>RectifiedDataset</wcs:CoverageSubtype>
            <wcs:nativeFormat>image/tiff</wcs:nativeFormat>
        </wcs:ServiceParameters>

    </wcs:CoverageDescription>

Each tag block represents a single image. The images are sorted by sensing time in descending order, i.e. the newest image is on top of the file. So we look at the first wcs:CoverageDescription and extract the tag wcs:CoverageId. In this case it is.

<wcs:CoverageId>RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES</wcs:CoverageId>

Note:

• The format of CoverageId is subject to change within the beta phase.
• Since the catalog is constantly updated, it is likely to find different and newer coverage ids.

To download the image from the server, we need the request method getCoverage and this coverage id.

curl -H "Authorization: Bearer ACCESS_TOKEN" -o /tmp/Antwerp_image_1.tiff -X GET "https://api.yaas.io/sap/earth-analysis/v1/wcs/getCoverage?coverageId=RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161124T175514_R008_V20161124T104352_20161124T104352T31UES"

If the image is not yet available on the server the request automatically triggers a background process, which retrieves the necessary data and calculates the requested image.
Then we get a short XML exception and HTTP return code 202; retry after a minute.

But what is this? The image is cloudy and incomplete!

Indeed clouds are a big issue in earth observation such as the multi-spectral images sensed by Sentinel-2 satellite. There is nothing we can do about it, except for better filter options. The downloaded XML query result contains information about cloud coverage in images. The tag cloudPercentage contains the percentage. The higher the value, the more cloudy is the image.

<om:result>
    <cloudPercentage uom="%">100.0</cloudPercentage>
</om:result>

Furthermore the images to download are only a part of a larger image. This larger images is sliced into fixed bounding boxes (each 100km x 100km), called tiles. Tiles on the border of the larger image result in images which are not filled completely.

To query only images which are less cloudy, we can extend the former query by dimensionTrim=cloudPercentage. We want to see only images which are maxmial 40 percent cloudy.

curl -H "Authorization: Bearer ACCESS_TOKEN" -o /tmp/Antwerp_query_result.xml -X GET "https://api.yaas.io/sap/earth-analysis/v1/wcs/describeEoCoverageSet?eoid=DS_RGB&dimensionTrim=lat(51.2,51.4)&dimensionTrim=long(4.1,4.6)&dimensionTrim=cloudPercentage(0,40)

Now the queried result contains less cloudy images. We will try an images from Sep 28th 2016 with coverage id RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161001T021620_R051_V20160928T105022_20160928T105637T31UES.

Lets try this one:

curl -H "Authorization: Bearer ACCESS_TOKEN" -o /tmp/Antwerp_image_2.tiff -X GET "https://api.yaas.io/sap/earth-analysis/v1/wcs/getCoverage?coverageId=RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161001T021620_R051_V20160928T105022_20160928T105637T31UES"

Note:
For this blog post, this image was reduces in size and quality. The orginal image is uncompressed and has a resolution of 10980×10980 pixel, where a single pixel represents 10m x 10m area. To get a better impression of the quality of Sentinel-2 images the area around the port is clipped (but still with lossy compression).

Where is what on the image?

Now having the image, it is important to know where the objects are located. Therefore we will transform pixel coordinates into world coordinates.

For those who have not heard about Spatial Reference Systems (SRS), yet, I recommend this article.

Since the image itself is planar and in a cartesian coordinate system, it is quite reasonable to use a planar cartesian coordinate system to reference pixel to world coordinates. The downloaded image is in UTM Zone 31U coordinate system, which refer to the EPSG code 32631.

Other images might refer to a different EPSG code.

Latter code can be found in the downloaded XML from the original query. Look for the tag gml:origin.

<gml:origin>
    <gml:Point srsName="http://www.opengis.net/def/crs/EPSG/0/32631" gml:id="RGB_________ESAS2___S2A_OPER_PRD_MSIL1C_PDMC_20161001T021620_R051_V20160928T105022_20160928T105637T31UES_grid_origin">
    <gml:pos>499980 5700000</gml:pos></gml:Point>
</gml:origin>
<gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/32631">10 0</gml:offsetVector>
<gml:offsetVector srsName="http://www.opengis.net/def/crs/EPSG/0/32631">0 -10</gml:offsetVector>

From this section we retrieve:

• EPSG code of SRS.
• Upper left corner coordinates of image.
• Size/dimension of pixels.

In short, for each pixel to the right, the position on earth moves 10 meters to east, and for each pixel to the bottom, the position moves 10 meters south. This SRS has further a unit of measure in meters, which means moving the position by 10 meters will change the coordinates by 10 units.

This said, the image pixel X=100, Y=250 (relative to upper left corner) means 1000m east and 2500m south. Relative to upper coordinate, this results in UTM Zone 31U coordinates 500980, 5697500.

But how to get the latitude and longitude from these values?

Therefore a projection library such as proj4 can be used. The respective python pachage is called pyproj and can be easily installed via pip.

We know the source SRS (EPSG 32631) and use a target SRS, whose EPSG code is 4326 (a widely used round earth SRS to use latitude and longitude).

import pyproj

SRS_TARGET = pyproj.Proj("+init=epsg:4326")

def transform_pixel_coordinate_to_long_lat(
        x,
        y,
        x_size=10,
        y_size=-10,
        upper_left_x=499980,
        upper_left_y=5700000,
        source_epsg=32631):
    source_srs_coordinate = (
        upper_left_x + x * x_size,
        upper_left_y + y * y_size)
    source_srs = pyproj.Proj(
        "+init=epsg:%i" % source_epsg)
    return pyproj.transform(
        source_srs,
        SRS_TARGET,
        source_srs_coordinate[0],
        source_srs_coordinate[1])

print transform_pixel_coordinate_to_long_lat(100, 250)
# (3.014096461455405, 51.42870180700858)

Lets check the coordinate transformation!

We use an image tool like Gimp to get pixel coordinates. We select a significant area for instance a corner in the northern part of the port (marked with the circle). Its pixel coordinates are X=08871 and Y=1616.

We use the python script to get the word coordinates of this pixel:

print transform_pixel_coordinate_to_long_lat(8871, 1616)
 # (4.272153046282123, 51.29896283926485)

Now swap the tuple and enter 51.29896283926485, 4.2721530462821231
in  Open Street Map and see where the pixel is located-

Congratulations! We just successfully downloaded an impressive image and found out world coordinates of an interesting pixel.

Summary

• We learned how to find out what bounding box we need and how to use it in the describeEoCoverageSet request method in images queries.
• Additionally we also filtered for less cloudy images.
• We downloaded images with the getCoverage request method.
• In the end we mapped pixel coordinates to world coordinates  (latitude and longitude).

With your new perspective on earth, you don’t need a window seat in the airplane anymore. So get started and write cool software.

Updates

• Dec 19th 2016: Renamed dimensionTrim parameter and respecitive XML tag to “cloudPercentage“