[GIS] Extracting data from a raster

gdalpythonrasterrasteriorioxarray

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

I'm not familiar with gdal and rasters and I'm using Python to extract the data out of the raster. Can you provide me the most efficient way to extract the data?

  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)

I get some value from the raster but the value is not correct.

My guess is that I am not converting the coordinates in the right way. I need to convert the coordinates 43.169129 lat, 15.174866 lon to coordinates used in the map to extract data at that exact spot.

Best Answer

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]
Related Question