Extracting data from a raster

Question

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

2. 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.

Answer 1

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]

Answer 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)

Answer 3

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