# 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

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)
AREA_OR_POINT=Area
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
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))

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.

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

## Related Questions

### Using FeatureCollection in Configurable Map Viewer?

1  Asked on August 18, 2021 by geebengs

### AttributeError: ‘Layer’ object has no attribute ‘metadata’ from ArcPy with ArcGIS Pro

2  Asked on August 18, 2021

### How to find Nearest point on line from a point cordinate in SQL server

1  Asked on August 18, 2021 by sanjith

### QGIS – How to render QgsMapCanvas on QApplication’s main window

1  Asked on August 18, 2021 by tilt_y

### Extracting Vertices / Snapping to points from Vector Tile Service (VTS) Layer

1  Asked on August 18, 2021

### gdalsrsinfo conflict between command line and python console

1  Asked on August 18, 2021

### Where can I get British Agricultural Field Boundaries as Vector Data?

2  Asked on August 18, 2021 by mathew-bayley

### GDAL code does not work on anaconda interpreter

2  Asked on August 18, 2021

### Shapefile PRJ to PostGIS SRID lookup table?

6  Asked on August 18, 2021 by ryankdalton

### How to calculate the accuracy of georeferenced aerial images in meters if the mean error is in pixels?

0  Asked on August 18, 2021 by pepe-pesh

### Selecting basins with discharge direction to East using ArcMap 10.3

0  Asked on August 18, 2021 by abdullah-nagy

### Qgis Work-around geometry-loss-mix-up-of-a-spatialite-view-joining-two-tables

0  Asked on August 17, 2021 by guest-no-1534

### Print layout of multiple maps revert to one identical map after save as image/zoom

0  Asked on August 17, 2021 by susan-miller

### Extract population count from GeoTIFF based on mask layer

2  Asked on August 17, 2021

### Mappa.js Clear Canvas with Leaflet marker

0  Asked on August 17, 2021

### Build tile server manually with PostgreSQL db on a different server

1  Asked on August 17, 2021 by vaishnavi-mukundhan

### Viewing TPW HDF4 (HDF-EOS4) Files in ArcGIS

0  Asked on August 17, 2021

### QGIS 3.10 – Changing display name of layer

1  Asked on August 17, 2021

### Maximum GPS-based location Sampling rate of smartphones

1  Asked on August 17, 2021 by jan-jan

### Building terrian data in ArcGIS Server for Cesium

0  Asked on August 17, 2021