2

I'm working with the ESA CCI land cover rasters. I'm trying to window and project them, but the results aren't what I expect.

If I run ncdump -h -s on one of the ESA CCI raster NetCDFs I get:

crs:wkt = "GEOGCS[\"WGS 84\", \n",
        "  DATUM[\"World Geodetic System 1984\", \n",
        "    SPHEROID[\"WGS 84\", 6378137.0, 298.257223563, \
             AUTHORITY[\"EPSG\",\"7030\"]], \n",
        "    AUTHORITY[\"EPSG\",\"6326\"]], \n",
        "  PRIMEM[\"Greenwich\", 0.0, AUTHORITY[\"EPSG\",\"8901\"]], \n",
        "  UNIT[\"degree\", 0.017453292519943295], \n",
        "  AXIS[\"Geodetic longitude\", EAST], \n",
        "  AXIS[\"Geodetic latitude\", NORTH], \n",
        "  AUTHORITY[\"EPSG\",\"4326\"]]" ;

I know also that the raster is 'north-up'...in the metadata from that same ncdump it is written: geospatial_lat_units = "degrees_north".

Now that I know the raster bounds, orientation and coordinate reference system (WGS84 = EPSG:4326), I attempt to window and project the raster using GDAL. I wish to project the data into EPSG:3035 (= ETRS89 = Europe), so I look up the bbox extent on epsg.io/3035 and read the following:

Projected bounds:
1896628.62 1507846.05
4656644.57 6827128.02

WGS84 bounds:
-16.1 32.88
39.65 84.17

With this in mind, I go ahead and make an attempt with the following (Python) script:

import os

os.system('gdal_translate -of vrt -a_ullr -180 90 180 -90 -a_srs "EPSG:4326" HDF5:"ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.nc"://lccs_class esacci.vrt')
os.system('gdalwarp -overwrite -of vrt -te_srs "EPSG:4326" -te -16.1 32.8 39.65 84.17 esacci.vrt esacci2.vrt')
os.system('gdalwarp -overwrite -s_srs "EPSG:4326" -t_srs "EPSG:3035" esacci2.vrt esacci.tif')

A glance at the resulting TIFF makes it seem like everything is in order. Furthermore, when I run gdalinfo on esacci2.vrt the corner coordinates are as expected (by epsg.io/3035):

Coordinate System is:
GEOGCS["WGS 84" ….]

Corner Coordinates:
Upper Left  ( -16.1000000,  84.1700000) ( 16d 6' 0.00"W, 84d10'12.00"N)
Lower Left  ( -16.1000000,  32.8000000) ( 16d 6' 0.00"W, 32d48' 0.00"N)
Upper Right (  39.6500000,  84.1700000) ( 39d39' 0.00"E, 84d10'12.00"N)
Lower Right (  39.6500000,  32.8000000) ( 39d39' 0.00"E, 32d48' 0.00"N)
Center      (  11.7750000,  58.4850000) ( 11d46'30.00"E, 58d29' 6.00"N)

HOWEVER, when I run gdalinfo on esacci.tif...the corner coordinates do not match what is predicted by epsg.io/3035:

Coordinate System is:
PROJCS["ETRS89 / LAEA Europe" ….]

Corner Coordinates:
Upper Left  ( 1894132.863, 6827128.020) ( 71d33' 9.14"W, 68d52'45.48"N)
Lower Left  ( 1894132.863, 1086761.052) ( 14d56'32.27"W, 29d14'19.31"N)
Upper Right ( 7060793.737, 6827128.020) ( 93d51'27.16"E, 66d16' 0.52"N)
Lower Right ( 7060793.737, 1086761.052) ( 37d59'10.65"E, 28d17' 0.95"N)
Center      ( 4477463.300, 3956944.536) ( 12d41'37.11"E, 58d41'13.10"N)

I can further replicate the EPSG:4326 to EPSG:3035 transform described on epsg.io/3035 with OSR, where the following results in the correct/predicted coordinates...which again are not what esacci.tif lists:

from osgeo import osr

s_srs = osr.SpatialReference()
s_srs.ImportFromEPSG(4326)
t_srs = osr.SpatialReference()
t_srs.ImportFromEPSG(3035)

transform = osr.CreateCoordinateTransformation(s_srs, t_srs)

print(transform.TransformPoint(-16.1, 32.88))
print(transform.TransformPoint(39.65, 84.17))

What am I missing and/or doing wrong here?

Improve this question
2
  • Why VRT format? What happens if you use -of GTIFF? Why are you warping twice? you can get a box from the original with -projwin ULX ULY LRX LRY or with your -te, like on the 2nd warp, but supply the -t_srs as one operation.
    – Michael Stimson
    Commented Nov 23, 2017 at 0:32
  • Putting the gdalwarp statement all on one line results in correct projected corner coords (from gdalinfo)! But the resulting TIFF does not seem correct (half of EU is clipped, part of NA shows). Alternatively, using a gdal_translate -projwin instead of the first gdalwarp yields similar results (correct corner coords, incorrect extent after visual check)...very confused..it seems a tradeoff: correct physical extent or correct corner coords.
    – Manselpotamus
    Commented Nov 23, 2017 at 9:34

2 Answers 2

Reset to default
3

You simply miss an 8 in -te -16.1 32.8 39.65 84.17: according to the bounds you have reported, it should be 32.88, not 32.8 in the first gdalwarp command.

Although, there's a more elegant way to convert & transform your data using gdalwarp as library function:

from osgeo import gdal

ds_in = gdal.Open('HDF5:"ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.nc"://lccs_class')
ds_out = gdal.Warp('esacci.tif', ds_in, 
                   format = 'GTiff',
                   dstSRS = 'EPSG:3035',
                   outputBounds = [-16.1, 32.88, 39.65, 84.17],
                   outputBoundsSRS = 'EPSG:4326')

# if you want to check the result
print gdal.Info(ds_out)

ds_out = None
Improve this answer
1
  • Thanks for this, it is true I forgot an 8! But that doesn't explain the magnitude of discrepancy I reported. In the end the issue was conceptual/trivial/basic: projected bounds are not the same thing as the corner coordinates of a rectangular raster..
    – Manselpotamus
    Commented Nov 27, 2017 at 11:32
2

The issue is conceptual and contradictory to the question the results were more or less as expected (minor discrepancy noted by Antonio Falciano's answer due to a missing hundred's place '8' in the bounds):

Projected bounds are not the same thing as the corner coordinates of a rectangular raster (after projection, if the result is not rectangular, missing value padding is added in order to create a rectangle..).

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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