# How to get value and matrix index at a given (x;y) point?
**URL:** https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692
**Category:** Geo
**Tags:** question, rasters
**Created:** [June 6, 2025, 9:24am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692 "2025-06-06T09:24:18Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [June 6, 2025, 9:24am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/1 "2025-06-06T09:24:18Z")
I have a set of rasters with the same grid and a large number of points and need to create a CZV with these points + information from the rasters (so my questions is similar to [this one](https://discourse.julialang.org/t/how-to-extract-raster-values-at-coordinate-locations/69752), but I am using the Rasters.jl package).
So I have a Raster object `myraster`:
```julia
julia> myraster
╭─────────────────────────╮
│ 669×609 Raster{Int16,2} │
├─────────────────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
↓ X Projected{Float64} LinRange{Float64}(2.985e6, 4.321e6, 669) ForwardOrdered Regular Points,
→ Y Projected{Float64} LinRange{Float64}(3.21e6, 1.994e6, 609) ReverseOrdered Regular Points
├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── metadata ┤
Metadata{Rasters.GDALsource} of Dict{String, Any} with 4 entries:
"units" => ""
"offset" => 0.0
"filepath" => "/home/lobianco/.julia/dev/GenFSM/dev_local/default/cache/res/fr/dtm/wc2.1_30s_elev_reg.tif"
"scale" => 1.0
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
extent: Extent(X = (2.985e6, 4.321e6), Y = (1.994e6, 3.21e6))
missingval: -32768
crs: PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
↓ → 3.21e6 3.208e6 3.206e6 3.204e6 3.202e6 3.2e6 3.198e6 3.196e6 3.194e6 … 2.01e6 2.008e6 2.006e6 2.004e6 2.002e6 2.0e6 1.998e6 1.996e6 1.994e6
2.985e6 -32768 -32768 -32768 -32768 -32768 -32768 -32768 -32768 -32768 258 264 268 272 279 276 270 287 344
⋮ ⋮ ⋱ ⋮
4.321e6 226 222 227 230 187 201 194 152 155 -32768 -32768 -32768 -32768 -32768 -32768 -32768 -32768 -32768
```
I have a point (x,y) using the same raster projections (3035):
```julia
julia> (x,y)
(3.9637574421965755e6, 2.949790197397309e6)
```
This looks inside my area (should be!), but if I try to retrieve the raster value with that values I have the following error:
```julia
julia> alt = myraster[X=At(x), Y=At(y)]
ERROR: SelectorError: attempt to select 3.9637574421965755e6 from lookup Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Float64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}} with bounds (2.985e6, 4.321e6)
```
Also, as I will have to lookup these values often, perhaps it is simpler to get the indexes and then just lookup the rasters with the matrix indexes.. how do I obtain the matrix indexes of a Raster object by providing the (x,y) coordinates using the same crs ?
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [June 6, 2025, 10:42am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/2 "2025-06-06T10:42:53Z")
I have found:
```julia
val_raster1 = myraster[Near(x), Near(y)]
```
or
```julia
Xdim = dims(myraster, Rasters.X)
Ydim = dims(myraster, Rasters.Y)
xcoords = collect(Xdim)
ycoords = collect(Ydim)
xidx = searchsortedlast(xcoords, x)
yidx = length(ycoords)- searchsortedlast(reverse(ycoords), y)
val_raster2 = myraster[xidx,yidx]
```
However they don’t give the same results.. by inspecting with QGIS, the second value is the correct one, while the first value (those with `Near()`) return the value of the eastbound cell… (actually I don’t know if this is a bug of Rasters.jl ??? )
EDIT:
Sometimes the `Near()` method returns the right pixel data, but it’s strange, it “fails” with coordinate points that are pretty at the center of their pixel..
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [June 6, 2025, 10:58am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/3 "2025-06-06T10:58:48Z")
`At` has to have exact values matching the lookup value, or within `atol` which you can pass with `At(val; atol=x)`. `Near(x)` should give you the closest point to `x`, west or east… I’m not sure what are doing with QGIS or if this is a bug, can you make a reproducible issue on github (file download/raster generation etc) for us to find out?
`Rasters.dims2indices(ras, (dims or lookups...,))` will give you the index if you need to use it multiple times.
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [June 6, 2025, 11:01am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/4 "2025-06-06T11:01:31Z")
But also note your Raster does not contain “pixels” that have a center (`Intervals` on x/y), but `Points`. So likely your confusion is from treating them as pixels starting from the point at one side.
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [June 6, 2025, 11:01am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/5 "2025-06-06T11:01:32Z")
Here a reproducible example:
```julia
import ArchGDAL
using Rasters, DimensionalData
tifraster = download("https://nc.beta-lorraine.fr/s/ga9Lexn6RJ9762g/download", "wc2.1_30s_elev_reg.tif")
myraster = Rasters.Raster(tifraster) # a 669x609 x→ ,y↓ raster
(x,y) = (4.214e6, 2.615e6)
val_raster1 = myraster[Near(x), Near(y)] # 1913
Xdim = dims(myraster, Rasters.X) # a DimensionalData.Sampled{<:Real} or similar
Ydim = dims(myraster, Rasters.Y)
xcoords = collect(Xdim)
ycoords = collect(Ydim)
xidx = searchsortedlast(xcoords, x) # 615
yidx = length(ycoords)- searchsortedlast(reverse(ycoords), y) #298
val_raster2 = myraster[xidx,yidx] # 1719 (QGIS value)
```
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [June 6, 2025, 11:02am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/6 "2025-06-06T11:02:46Z")
Thanks, would be great if you can paste to github, it will be forgotten here as its not in a formal queue
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [June 6, 2025, 11:04am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/7 "2025-06-06T11:04:11Z")
ok, I didn’t realize that.
I will post it on GitHub, but I now suspect the issue is with the data being indeed a point instead of (as I thought) a pixel…
Thank you for the very quick reply…
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [June 6, 2025, 11:06am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/8 "2025-06-06T11:06:19Z")
Yeah, points/intervals is a bit confusing - thats why its printed for every dimension.
Many GIS tools just ignore the point/pixel metadata and treat everything as pixels, but for `Near` thats not at all correct - the nearest point is not the same as the nearest pixel if that point is the left/top etc edge.
(you can use `intervalraster = set(x, Intervals)` to change your raster to intervals)
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [June 6, 2025, 11:10am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/9 "2025-06-06T11:10:08Z")
but in this case (as you can see from the screenshot) the coordinate to lookup is very close to the center of (what I thought was ) a pixel…
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [June 27, 2025, 5:04am UTC](https://discourse.julialang.org/t/how-to-get-value-and-matrix-index-at-a-given-x-y-point/129692/10 "2025-06-27T05:04:10Z")
Sorry I missed your reply. If you want it at the edge of the pixel, you need Intervals(Start()) as Center() is the default when converting from points.