Reprojections and great circles

GRASS GRASS user
, ,
1

I’m not sure if this is a GRASS issue or a QGIS issue, but I start here. Visually, there is no issue in GRASS, but there is an issue in QGIS.

I’m working with two grids and trying to reproject between the two and map cell areas via a third grid.

Grid E (red) is a global EPSG:4326 lat/lon grid that is built on the equirectangular projection, meaning the horizontal lines follow latitude bands and not the great circle.

Grid P (blue) is a polar grid on EPSG:3031.

I’m building an exchange grid X (black) which is the intersection of the two using v.proj and then v.overlay.

4326
4326871×688 201 KB

When I export to GeoPackage and view all three grids QGIS in EPSG:4326 everything looks correct. See image above, which shows all three grids (red E, blue P, and black X). The exchange grid (black X) is shown to perfectly match all the other grids here, lying directly underneath them.

3031
3031871×688 78.9 KB

When I view this same region in EPSG:3031, it does not look correct. The same large red E grid cell is shown above, and then the cell with the dot is in below.

3031_zoom
3031_zoom871×688 24 KB

Zoomed out things look OK, the black X grid is aligned under the other two. But zoomed in, it is clear that the X and E grids do not overlap. I believe the curved black X grid follows a latitude band, and the straight red E line is a great circle between the two corners of that grid cell.

If I view this same region in GRASS, all lines correctly overlap:

grass
grass979×781 51.6 KB

I don’t particularly care about the graphics. But I’m using GRASS “v.to.db map=X option=area units=meters columns=area” and I need the area to be computed correctly.

Does anyone know where the problem is here? Should this issue be raised with QGIS? With GDAL which I believe exports from GRASS to the GeoPackage? Somewhere else? Or no issue if area is correct, and this is a design decision of an underconstrained system that is intention and hopefully documented somewhere?

Thanks,

-k.

2

v.to.db computes the areas in latlon CRS using G_ellipsoid_polygon_area, in projected CRS with G_planimetric_polygon_area.

Somewhat related to that v.proj has smax parameter to more accurately reproject vectors by increasing the number of vertices.

Not sure about the graphics…

3

Thank you for those links and explanation.

From the source code for computing area,

This function correctly computes (within the limits of the series approximation) the area of a quadrilateral on the ellipsoid when two of its sides run along meridians and the other two sides run along parallels of latitude.

This makes me think this is not the right tool to use for vectors that are arbitrarily rotated or made of of complex shapes that do not have edges running along latitude or longitude lines.

I’ll do some comparisons between the GRASS area computation and Planimeter from GeographicLib · GitHub or some other tools.

4

See also (long discussion):

5

Thanks for this link. Helpful. Maybe another discord issue: I did not get an email about your 2nd reply, only your 1st. But I am “subscribed” to this topic.

6

Hmm, no idea, perhaps check email and notification settings. I don’t think I had that issue. Did you get this message?

7

Yes :). Thank you.