How to get value and matrix index at a given (x;y) point?

Specific Domains Geo
,
1

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, but I am using the Rasters.jl package).

So I have a Raster object myraster:

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> (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> 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 ?

2

I have found:

val_raster1 = myraster[Near(x), Near(y)]

or

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

3

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.

4

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.

5

Here a reproducible example:

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)

Screenshot from 2025-06-06 12-59-44
Screenshot from 2025-06-06 12-59-441109×732 12.2 KB

6

Thanks, would be great if you can paste to github, it will be forgotten here as its not in a formal queue

7

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…

1 Like
8

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)

9

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…