AnswerBun.com

Extracting data from a raster

Geographic Information Systems Asked by ktop on December 17, 2020

I’m working on a project where I must 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 get the data out of the raster. I’m so confused?

  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.

3 Answers

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]

Correct answer by snowman2 on December 17, 2020

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

Answered by Jose on December 17, 2020

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)

Answered by user97103 on December 17, 2020

Add your own answers!

Related Questions

Shapefile PRJ to PostGIS SRID lookup table?

6  Asked on August 18, 2021 by ryankdalton

     

Mappa.js Clear Canvas with Leaflet marker

0  Asked on August 17, 2021

 

Maximum GPS-based location Sampling rate of smartphones

1  Asked on August 17, 2021 by jan-jan

     

Ask a Question

Get help from others!

© 2022 AnswerBun.com. All rights reserved. Sites we Love: PCI Database, MenuIva, UKBizDB, Menu Kuliner, Sharing RPP