5

I'm currently working on a project where I have to use the Corine Land Cover (CLC) map. https://land.copernicus.eu/pan-european/corine-land-cover/clc2018 (You can download the map in the .tif format)

My problem is how to translate the real-world coordinates:

43.1691929 lat, 15.174866 lon

to coordinates used in the map to extract data at that spot. I got some value from the raster but the value is not correct.

  1. Gdal info:
Coordinate System is:
PROJCRS["ETRS_1989_LAEA",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101004,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["Lambert Azimuthal Equal Area",
        METHOD["Lambert Azimuthal Equal Area",
            ID["EPSG",9820]],
        PARAMETER["Latitude of natural origin",52,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",10,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",4321000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",3210000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",3035]]
Data axis to CRS axis mapping: 1,2
Origin = (900000.000000000000000,5500000.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  900000.000, 5500000.000) ( 56d30'18.51"W, 56d29' 4.75"N)
Lower Left  (  900000.000,  900000.000) ( 23d49'33.58"W, 24d17' 3.04"N)
Upper Right ( 7400000.000, 5500000.000) ( 72d54'22.09"E, 58d57' 9.90"N)
Lower Right ( 7400000.000,  900000.000) ( 40d39'45.75"E, 25d32'40.96"N)
Center      ( 4150000.000, 3200000.000) (  7d30'57.52"E, 51d53' 2.21"N)
Band 1 Block=65000x1 Type=Int16, ColorInterp=Gray
  Min=111.000 Max=999.000   Computed Min/Max=111.000,999.000
  Minimum=111.000, Maximum=999.000, Mean=326.518, StdDev=118.029
  NoData Value=-32768
  Metadata:
    DESCRIPTION=clc18
    RepresentationType=THEMATIC
    STATISTICS_MAXIMUM=999
    STATISTICS_MEAN=326.51842078382
    STATISTICS_MINIMUM=111
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=118.02878635921
    STATISTICS_VALID_PERCENT=24.58
  1. Python Code
import gdal
import numpy
from affine import Affine


lons=[15.174866]
lats=[43.169129]


fn="C:/path-to-the-map/map.tif"

ds=gdal.Open(fn)


transform=ds.GetGeoTransform()
xOrigin=transform[0]
yOrigin=transform[3]
pixelWidth=transform[1]
pixelHeight=transform[5]



aff=Affine.from_gdal(xOrigin,pixelWidth,0.0,yOrigin,0.0,pixelHeight)


x_coords,y_coords=aff*(numpy.array(lons),numpy.array(lats))

band=ds.GetRasterBand(1).ReadAsArray()


x=int(x_coords[0]/pixelWidth)
y=int(y_coords[0]/pixelHeight)


value=band[x][y]
print(value)

My guess is that I am not converting the coordinates in the right way.

1
  • Welcome to Geographic Information Systems! Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Your questions should as much as possible describe not just what you want to do, but precisely what you have tried and where you are stuck trying that. Please check out our short tour for more about how the site works.
    – Ian Turton
    Apr 12, 2020 at 16:20

3 Answers 3

15

You can use rioxarray for this:

import rioxarray
from pyproj import Transformer

# convert coordinate to raster projection
lon = 15.174866
lat = 43.169129

rds = rioxarray.open_rasterio("C:/path-to-the-map/map.tif")
transformer = Transformer.from_crs("EPSG:4326", rds.rio.crs, always_xy=True)
xx, yy = transformer.transform(lon, lat)

# get value from grid
value = rds.sel(x=xx, y=yy, method="nearest").values

You can also do this with rasterio:

import rasterio
from pyproj import Transformer

lon = 15.174866
lat = 43.169129

with rasterio.open("C:/path-to-the-map/map.tif") as rds:
    # convert coordinate to raster projection
    transformer = Transformer.from_crs("EPSG:4326", rds.crs, always_xy=True)
    xx, yy = transformer.transform(lon, lat)

    # get value from grid
    value = list(rds.sample([(xx, yy)]))[0]

0
2

If your expected output is just convert raster information to text, you can try gdal2xyz.py

$ python gdal2xyz.py -csv input.tif output.csv

Then you will get a csv file with XYZ information: Longitude, Latitude, and Raster information (land use)

0
2

Similarly, you can use gdallocationinfo from the command line. E.g. if you have your coordinates in a text file called coords.xy, you can just type

cat coords.xy | gdallocationinfo -wgs84 -valonly map.tif > values.xy

and your outputs are in file values.xy (same order as coords.xy). Blindingly fast! See here for details

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.