cellsize doesn't seem to be correct

[alex-s-gardner]

Here is an example of a north polar sterographic projected dataset with a map grid resolution of 120m. Scale distortion makes actual areas near the pole slightly larger and actual areas away from the poles slightly smaller than map projected areas. The distortion >65 deg latitude is < 5%... this would be the right hand side of the image that is shown in the example below. Map cell area is 120m*120m = 0.0144 km^2 but Rasters.cellsize returns areas > 0.028 or 2x map scale.

using Rasters
using Plots
path2rasterfile = "https://its-live-data.s3.amazonaws.com/velocity_mosaic/v2/static/cog/ITS_LIVE_velocity_120m_RGI01A_0000_v02_vx.tif"

foo = Raster(path2rasterfile)
ca = Rasters.cellsize(foo[1:10:end,1:10:end]) / 10^2
heatmap(ca)

[rafaqz]

Best not to use 10k * 20k files as an mwe ;)

Can you just post the plot? It might be because cellsize is on a perfect sphere?

[rafaqz]

Maybe @tiemvanderdeure knows whats going on here

[tiemvanderdeure]

cellsize is always a bit off because it assumes the world is spherical. It also reprojects everything to WGS84 in order to calculate the angle between the points, which I would assume is the problem here? But it seems weird that it should be off by this much, so I'm not entirely sure what's going on.

[tiemvanderdeure]

I doubt this would be because cellsize is so distorted at higher latitudes per se. I can calculate the area of greenland within 1% like this:

using NaturalEarth, Rasters
countries = naturalearth("admin_0_countries", 110)
greenland = countries[findfirst(contains("Greenland"), countries.NAME_EN)].geometry
greenland_ras = rasterize(greenland, res = 0.1, fill = true, missingval = false, crs = EPSG(4326))
greenland_size = sum(cellsize(greenland_ras)[greenland_ras])

Which makes me think that it must be the reprojection step that's the problem.

[alex-s-gardner]

for the MWE posted above this is the plot of Rasters.cellsize computed area vs cell area in map coordinates... the ratios is the area scale distortion which should not be < 0.9 as one move toward the pole.. right hand side of image

map_cellarea_km2 = (120*120)/(1000 .^ 2)
heatmap(map_cellarea_km2 ./ ca)

[alex-s-gardner]

Here is a plot of the EPSG 3413 scale distortion in x and y as a function of latitude (here 1m in x/y on ground is the distance on ground assuming the WGS84 ellipsoid)

[alex-s-gardner]

Yup... projected area is off by a factor of 3

greenland_size for EPSG 4326 = 2.2E6 greenland_size for EPSG 3413 = 6.1E6

Also very slow for EPSG 3413.

using NaturalEarth, Rasters, GeometryOps
countries = naturalearth("admin_0_countries", 110);
greenland = countries[findfirst(contains("Greenland"), countries.NAME_EN)].geometry;
greenland_ras = rasterize(greenland, res=0.1, fill=true, missingval=false, crs=EPSG(4326));
greenland_size_4326 = sum(cellsize(greenland_ras)[greenland_ras])

greenland = GeometryOps.reproject(greenland; source_crs=GFT.EPSG(4326), target_crs=GFT.EPSG(3413), calc_extent=true);
greenland_ras = rasterize(greenland, res=500, fill=true, missingval=false, crs=EPSG(3413));
greenland_size_3413 = sum(cellsize(greenland_ras)[greenland_ras])

[tiemvanderdeure]

Okay, I think I figured it out. Cellsize reprojects via GDAL, and GeometryOps and ArchGDAL disagree about the order the X and Y coordinates should be in.

julia> GeometryOps.reproject([-3e6, 0]; source_crs = EPSG(3413), target_crs = EPSG(4326))
(-135.0, 62.80769951036124)

julia> ArchGDAL.reproject([-3e6, 0], EPSG(3413), EPSG(4326))
2-element Vector{Float64}:
   62.80769951036124
 -135.0

I'm not sure which projections cellsize is broken for. We do have a test against some equal area projections.

Maybe we should switch to using GeometryOps/Proj?

[rafaqz]

Oh you just have to pass a keyword to ArchGDAL to fix it to traditional order.

Probably Rasters does that?

[tiemvanderdeure]

It seems like that keyword changes the order of the input point though, the transformed point is still lat-lon (Y-X) for EPSG4326.

Never mind, you were right. We should use that.