Function [Latitude & Longitude] -> Index

[AquaPore]

Hi,

I was looking for a function in Rasters which converts Longitude and Latitude into index [iX, iY] of a map.tiff.

As my knwoledge in Rasters is limited I wrote one see below.

I am writing to determine if there is a function available in Rasters.jl which does the same thing so that I can compare the results.

Many thanks for any suggestions,
Joseph

uisng Rasters

Map = Rasters.Raster("Map.tiff")

function LAT_LONG_2_INDEX(;Map, [146705.520 , 42165.261])
         Longitude_X = Param_GaugeCoordinate[1]
         Latitude_Y  = Param_GaugeCoordinate[2]

         Longitude   = Rasters.lookup(Map, X)
         Latitude    = Rasters.lookup(Map, Y)

         Nlongitude  = length(Longitude)
         Nlatitude   = length(Latitude)

			# Longitude
				ΔX = Longitude[2] - Longitude[1]
				@assert(Longitude[1] - ΔX / 2.0 ≤ Longitude_X ≤ Longitude[Nlongitude] + ΔX / 2.0)
				iX = 1
				for i=1:Nlongitude
					if Longitude[i] - ΔX / 2.0 ≤ Longitude_X < Longitude[i] + ΔX / 2.0
						iX = i
						break
					end
				end # for i=1:Nlongitude

			# Latitude
			 	ΔY = Latitude[2] - Latitude[1]
				@assert( Latitude[Nlatitude] + ΔY / 2.0  ≤ Latitude_Y ≤ Latitude[1] - ΔY / 2.0)
				iY = Nlatitude
				for i=Nlatitude:-1:1
					if Latitude[i] + ΔY / 2.0  ≤ Latitude_Y < Latitude[i] - ΔY / 2.0
						iY = i
						break
					end # if Latitude_Sort[iY] ≥ Lat_Y
				end # i=1:Nlatitude

				println( "LAT_LONG_2_INDEX: [$iX ; $iY]" )
		return iX, iY
		end

[asinghvi17]

DimIndices(ras)[X(At(lon)), Y(At(lat))] should give you a tuple (X(iX), Y(iY)), you can call val on those elements to get the integers. And you can use any selector.

But you can also just do ras[X(At(lon)), Y(At(lat))] to get the value at those indices.

[AquaPore]

Thanks @asinghvi17 for your kindness.

I finally understood how At() works.

QUESTION 1

Your proposed method works very well for all the Vector of lon and lat. But for some reason does not seem to work for a single point

ras[X(At(146704.657)), Y(At(42164.757))]

I get the following error

SelectorError: attempt to select 146704.657 from lookup
 Projected{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, 

May be one should use
ras[X(Near(146704.657)), Y(Near(42164.757))]

QUESTION 2

I would like to get a function for example which returns the index [iX, iY ] of the ras at the coordinate [146704.657, 42164.757], how do I do?

Many thanks for your kind help,
Joseph

[asinghvi17]

Ah yeah, Near is better for floating point.

Did you try the DimIndices method I mentioned? That should get you the linear indices. Otherwise there is also dims2indices(ras, (X(...), Y(...)))

[AquaPore]

Thanks @asinghvi17,

You are ever so kind and knowledgeable.

It works and the solution is close to the result of my code +/- 1 which I will investigate.

iX, iY = Rasters.dims2indices(Map, (X(Rasters.Near(146704.657)), Y(Rasters.Near(42164.757))) )
 [361, 47]

May be it will be usefull for others to update the manual.

Question

Would you know how to do the oposite to double check:

Function[iX, iY] = Coordinate [146704.657, 42164.757] with an accuracy of cell size.

With kindness,
Joseph

[asinghvi17]

DimPoints(ras)[iX, iY] ought to do it! Assuming x is your first dimension and y is your second. You can safeguard against this not being the case by

DimPoints(ras)[X(iX), Y(iY)]

P.S. That ±1 could be because of sampling (Points vs Intervals{Start} vs Intervals{End}) which may not have been accounted for in the original code. You could try shiftlocus to switch the intervals around from start to end or vice versa, there should be some documentation on that in DimensionalData

[AquaPore]

Thank you for your kindness,

This worked:
DimPoints(Gauge)[X(iX), Y(iY)]

I will check now with QGIS by creating a point in a raster format and I will test by using different resolutions. But it seems tricky to solve the +/- 1 issue, which is documented in
how-to-get-value-and-matrix-index-at-a-given-x-y-point

I will keep you posted,

With kindness,
Joseph

[AquaPore]

Dear Anshul Singhvi,

A great thanks for your kindness, for my education, I wrote LAT_LONG_2_INDEX() and I am happy that it gave same results as your solution.

iX, iY = Rasters.dims2indices(Map, (X(Rasters.Near(146704.657)), Y(Rasters.Near(42164.757))) )
 [361, 47]

As suspected I raised a BUG: Point coordinate -> Raster point is shifted
#997

I am now closing the issue,

With kindness,
Joseph Pollacco

function LAT_LONG_2_INDEX(;Map, Param_GaugeCoordinate)
         Longitude_X = Param_GaugeCoordinate[1]
         Latitude_Y  = Param_GaugeCoordinate[2]

         Longitude   = Rasters.lookup(Map, X)
         Latitude    = Rasters.lookup(Map, Y)

         Nlongitude  = length(Longitude)
         Nlatitude   = length(Latitude)

			# Longitude
				ΔX = Longitude[2] - Longitude[1]

				@assert(Longitude[1] ≤ Longitude_X ≤ Longitude[Nlongitude])
				iX = Nlongitude
				for i=1:Nlongitude
					if Longitude[i] - ΔX / 2.0 ≤ Longitude_X < Longitude[i] + ΔX / 2.0
					# if Longitude[i] - ΔX ≤ Longitude_X < Longitude[i]
						iX = i
						break
					end
				end # for i=1:Nlongitude

			# Latitude
			 	ΔY = Latitude[2] - Latitude[1]
				@assert( Latitude[Nlatitude] + ΔY / 2.0  ≤ Latitude_Y ≤ Latitude[1] - ΔY / 2.0)
				iY = 1
				for i=Nlatitude:-1:1
				if Latitude[i] + ΔY / 2.0  ≤ Latitude_Y < Latitude[i] - ΔY / 2.0
					# if Latitude[i]  < Latitude_Y ≤ Latitude[i] - ΔY
						iY = i
						break
					end # if Latitude_Sort[iY] ≥ Lat_Y
				end # i=1:Nlatitude

				println( "LAT_LONG_2_INDEX: [$iX ; $iY]" )
		return iX, iY
		end

[AquaPore]

A great thanks to @asinghvi17,

This was really helpfull, and it will be usefull to update the manual to put to light the amazing functionality.

iX, iY = Rasters.dims2indices(Map, (X(Rasters.Near(146704.657)), Y(Rasters.Near(42164.757))) )

[lazarusA]

well all of these as selectors, PR-suggestions on improving the corresponding section are always welcome! https://rafaqz.github.io/Rasters.jl/dev/get_started#Subsetting-an-object

See also: https://rafaqz.github.io/DimensionalData.jl/stable/selectors for more examples.