BUG: Point coordinate -> Raster point is shifted

[AquaPore]

########################

WHAT WE WANT TO ACHIEVE

In order to delineate the catchment area we need to accurately provide the location coordinate of the gauging station e.g. [146709.6323,42159.6226]. We also need to assure that the location of the gauging station is exactly located on a rasterized river network.

I would have imagined that the location of the gauging station should be located at the centre of the cell. Nevertheless, this method was giving me random results which shifted the cell Up, Bottom, Left, Right.

I investigated the problem and found that depending on the point location of the gaging station in the cell (displayed in QGIS) would give different results (see diagram):
• Left-top;
• Right-top;
• Left-Bottom;
• Right-Bottom;

In order that the location of the gaging station matches the location of the raster cell, I found that the location of the gaging station should be at the left-top of the cell (In the Northern hemispher) and not in the centre of the cell, which I believe is a BUG.

########################

# CONVERTING GAGING POINT -> RASTER POINT
function GAUGE(;Gauge_Name, Longitude=Longitude, Latitude=Latitude, Param_GaugeCoordinate)

	Gauge = Rasters.Raster((Longitude, Latitude), crs=Crs_GeoFormat)
	Gauge .= 0

		iX_Gauge, iY_Gauge = Rasters.dims2indices(Gauge, (X(Rasters.Near(Param_GaugeCoordinate[1])), Y(Rasters.Near(Param_GaugeCoordinate[2]))))

		Gauge[iX_Gauge, iY_Gauge] = 1

	# Saving file
		Path_OutputGauge = joinpath(Path_Root, Path_OutputWflow, Gauge_Name)
		Rasters.write(Path_OutputGauge, Gauge; ext=".tiff", force=true, verbose=true, missingval=0)
end

# PARAMETERS
ΔX       = step(dims(Map, X))
ΔY       = step(dims(Map, Y)

Coord_X_Left   = first(dims(Map, X))
Coord_X_Right  = last(dims(Map, X))
Coord_Y_Top    = first(dims(Map ,Y))
Coord_Y_Bottom = last(dims(Map,Y))
Crs_GeoFormat = GeoFormatTypes.convert(WellKnownText, EPSG(Param_Crs))

Longitude, Latitude = Rasters.X(Coord_X_Left: ΔX: Coord_X_Right, crs=Crs_GeoFormat), Rasters.Y(Coord_Y_Top: ΔY: Coord_Y_Bottom, crs=Crs_GeoFormat)

# RUNNING CODE

Param_GaugeCoordinate_X_Right_Top    = [146709.6323,42159.6226]
Param_GaugeCoordinate_X_Left_Top     = [146700.2167,42159.7300]
Param_GaugeCoordinate_Center         = [146726.378,42150.827]
Param_GaugeCoordinate_X_Right_Bottom = [146709.7495,42150.1972]
Param_GaugeCoordinate_X_Left_Bottom  = [146700.2801,42150.2118]

GAUGE(Gauge_Name="Gauge_Wflow_X_Right_Top.tiff", Param_GaugeCoordinate=Param_GaugeCoordinate_X_Right_Top)
GAUGE(Gauge_Name="Gauge_Wflow_X_Left_Top.tiff", Param_GaugeCoordinate=Param_GaugeCoordinate_X_Left_Top)
GAUGE(Gauge_Name="Gauge_Wflow_X_Right_Bottom.tiff", Param_GaugeCoordinate=Param_GaugeCoordinate_X_Right_Bottom)
GAUGE(Gauge_Name="Gauge_Wflow_X_Left_Bottom.tiff", Param_GaugeCoordinate=Param_GaugeCoordinate_X_Left_Bottom)

[rafaqz]

This isn't a bug, it's simply different across different backends (GDAL vs CF standards) where the point is given within a pixel. Corner or center.

When defining it yourself, you can specify you want intervals sampling with the "locus" in the corner. This is X(vals; sampling=Intervals(Start())) and the same on Y. Center() means the point is the center of the pixel, End() (rarely used) is the other corner.

(Also, by default you get Points rather than Intervals which by definition here have Center locus, and they keep that when converted to Intervals. There just needs to be some default and for better or worse that's what it is. In CF world even dates are always at interval centres, as painful as that is to use, so it's not without reason)

[AquaPore]

Dear @rafaqz,

A great thanks for your prompt reply, I have now a better understanding on the challange.

So if I got in roight, in Julia to implement the different sampling method:

	Gauge = Rasters.Raster((Longitude, Latitude), crs=Crs_GeoFormat)
	Gauge .= 0

        # Locus in center
		Gauge = Rasters.set(Gauge, Rasters.Center)
          # Corner left
		Gauge = Rasters.set(Gauge, X => Rasters.Intervals(Rasters.Start()))
          # Corner right
		Gauge = Rasters.set(Gauge, X => Rasters.Intervals(Rasters.End()))

I am expecting to get the same result [iX_Gauge, iY_Gauge] when using Rasters.set(Gauge, Rasters.Center) independently on the accuracy of the point we select GaugeCoordinate_X_Right_Top, GaugeCoordinate_X_Left_Top, GaugeCoordinate_X_Right_Bottom, ' The cell size is 10 x 10m

		iX_Gauge, iY_Gauge = Rasters.dims2indices(Gauge, (X(Rasters.Near(Param_GaugeCoordinate[1])), Y(Rasters.Near(Param_GaugeCoordinate[2]))))

I will be ever greatfull, if you can let me know if I missed something?

With kindness,
Joseph

[rafaqz]

I don't totally understand your example, but you have only set the X axis? So the Y will still be Center for all those examples. The axes are totally independent.

With Intervals you should also prefer to use Contains rather than Near as it makes sure your point is actually inside an interval.

[AquaPore]

Dear @rafaqz ,

A great tthanks for all your help you provided, and appologize if I was not so clear.

I am questioning what will be your recommended practice to specify the interval sampling method (e.g. Center(), End(), Start()) at the very beging when for example we are reading a map:

Dem = Rasters.Raster(Path; Center())

If this is not possible, I am still confused what will be your recommended method of setting the interval sampling method of a map?

With kindness,
Joseph

[rafaqz]

Normally when you load a file the locus should be correct?

Before trying random ideas lets get clear what is happening:

• Do you think the file was written incorrectly?
• Was it a gdal tiff or netcdf/zarr/grid file?
• Does it load as Points or Intervals? (feel free to paste the printout for the raster here)

[rafaqz]

Also, a note on making this easier for me to understand your problem I'm struggling a bit and havent really understood the problem yet:

• when you make an example always try printing the results so we can see what the raster settings are
• if you include file downloads so your script just runs I will just run it and tell you the problem very quickly. If it doesn't run it may take months!
• you can code however you like but I honestly struggle to read non-standard julia formatting! like caps in names are usually types, indenting needs to be consistent and lines shouldn't go off the screen in github. try https://github.com/JuliaDiff/BlueStyle

[rafaqz]

Oh actually I think your problem is just using Near rather than Contains! Near is perfect for Points but is unreliable for Intervals as it just the closest to the point in the lookup value. Maybe it shouldn't be... But if you have intervals what you mean is Contains, i.e. you want the interval that contains that value.

[AquaPore]

Dear @rafaqz,

A great thanks for pointing me to the right direction by recommending using Near rather than Contains.

I had a closer look at your excellent manual Sinusoidal-Projection-MODIS and I think I found what I was looking for:

using DimensionalData, Rasters
Dem     = Rasters.Raster(Path)
Dem_Resample = DimensionalData.shiftlocus(Rasters.Center(), Dem)

A geat thanks for your kind help,
Joseph

[AquaPore]

Dear @rafaqz,

May be this is out of topics of Julia style. I have a long code, and I personally find easy to read when the functions I create are all in caps. But I undestand that this style is not universally accepted.

function DEM_DERIVE_COASTLINES(;Dem, Longitude, Latitude, Crs, Missing=NaN, ϵ=0.001, DemMin=0.5)

I would try and provide you the input file when I raise an issue.

Kind regards,
Joseph

[rafaqz]

Yeah - there are a million good ways to write code! But its about what we all see every day being consistent.

I personally struggle to read code that isn't formatted like all the other code I see. When there's a very short time available to understand an issue this really matters - it will take me a few times longer to understand and often that means I just wont ever understand it in the time I can allocate to it.