New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
PROJ6+ + WKT2 support #1240
PROJ6+ + WKT2 support #1240
Conversation
This is huge! Tested with importing vectors in EPSG:3059 without reprojecting and also with reprojecting on import into EPSG:4326 – works fine. Didn't test places where improper axis ordering might pop-up. At least in importing and displaying everything is fine. |
Huge work @metzm ! The CI shows two errors with Ubuntu 16 (old but maybe still fixable):
and
Oh, it is an antique GDAL: |
+1 for dropping Ubuntu 16 |
Note that also ubuntu-18.04 fails. There is a compiler macro condition block in
which causes failure with GDAL < 3. |
Amazing work! So far tested the first suggested test with PROJ 7.2.1 and GDAL 3.2.1. Seem to work just fine! Compiling with clang issued the following three warnings on files you modified. They are not caused by this PR, but maybe you wouldn't mind addressing them while you're at it. lib/proj/do_proj.c :
raster/r.in.gdal/main.c :
vector/v.out.ogr/main.c :
|
I don't think that would be necessary, at least not for this enhancement. But this needs a good testing with GDAL <3 |
I encountered problems with projecting back and forth between epsg 23700 and 3035. #!/usr/bin/env bash
GISDBASE=$HOME/grassdata
GRASS=`which grass79`
cat << EOF > script.sh
#!/usr/bin/env bash
v.proj location=projtest_3035 mapset=PERMANENT input=point_23700_2 output=point_23700_2
echo "original:"
v.out.ascii input=point_23700@PERMANENT format=point
echo "back-and-forth:"
v.out.ascii input=point_23700_2@PERMANENT format=point
EOF
chmod +x run.sh
cat << EOF > point_23700.geojson
{
"type": "FeatureCollection",
"name": "point_23700",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::23700" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 737000.0, 240000.0 ] } }
]
}
EOF
$GRASS -c --text EPSG:23700:3 $GISDBASE/projtest_23700 --exec v.import input=point_23700.geojson
$GRASS -c EPSG:3035 $GISDBASE/projtest_3035 --exec v.proj location=projtest_23700 mapset=PERMANENT input=point_23700 output=point_23700_2
$GRASS --text $GISDBASE/projtest_23700/PERMANENT --exec $(pwd)/script.sh results in:
and
While using ogr2ogr is pretty straightforward:
|
Fine with me. Fair enough for 8.0 and a major improvement like this one. |
Add missing parenthesis Co-authored-by: nilason <n_larsson@yahoo.com>
* fix some compiler warnings * fix support for GDAL < 3 * fix vector reprojection of a single point
Thanks all for your helpful comments! I have fixed support for older GDAL/PROJ versions. Removing support for older GDAL/PROJ versions would be very helpful! The code base of the respective libraries and modules could be substantially cleaned up (removing lots of these Compiler warnings have been taken care of, and reprojection of a single point works now. As long as my original test fails, this PR will remain a draft. I am looking into it. |
Ok, see #1254 |
The last commit ab7a768 fixes issues with the area of interest to be set for PROJ 6+. However, my tests still fail. |
There is a bug in PROJ causing my tests to fail with |
This PR (together with previous commits) would fix bugs also present in G78: CRS mismatch between input and output when creating a location from some GDAL/OGR recognized input, importing to this location and exporting again. The CRS should be exactly preserved. This bug happens when GRASS is compiled with PROJ6+ and GDAL3. A bug in PROJ can cause PROJ to select a wrong coordinate operation with a datum transformation grid that does not cover the whole area of interest. In these cases, reprojection fails. This PR takes care of this bug. Fixing these bugs in G78 requires that all changes related to the new SRID and WKT2 handling including all affected import and export modules need to be backported. |
My previous observations, issues have all been solved.
This may be an issue for a follow-up PR though. |
This looks like a GUI issue. Please provide the full GUI output. |
Not reproducible through the menu "Raster -> Develop raster map -> Reproject raster map from different location". Must be a bug in the Data tab. |
Thanks very much @rouault for your suggestions! |
Yes, its only in Data tab and also unrelated to this PR. Reported with #1273. |
* free memory * code cleanup
No possibility for me to test this seriously in the near future, but big kudos to @metzm for this great improvement ! |
Co-authored-by: nilason <n_larsson@yahoo.com> * adjust lib/gis, lib/proj, g.proj, r.proj, v.proj, r.in.gdal, r.external, r.out.gdal, v.in.ogr, v.external, v.out.ogr
@neteler backported to G78 |
Co-authored-by: nilason <n_larsson@yahoo.com> * adjust lib/gis, lib/proj, g.proj, r.proj, v.proj, r.in.gdal, r.external, r.out.gdal, v.in.ogr, v.external, v.out.ogr
* r.in.nasadem + r.in.srtm.region: fix PROJ parsing Due to OSGeo/grass#1240 the output structure of `g.proj -g` changed, leading to ``` r.in.nasadem user='neteler' password='XXXXXXXXX' output=nasadem resolution=30 Traceback (most recent call last): File "/home/mundialis/.grass7/addons/scripts/r.in.nasadem", line 631, in <module> main() File "/home/mundialis/.grass7/addons/scripts/r.in.nasadem", line 475, in main SRCGISRC, TMPLOC = createTMPlocation() File "/home/mundialis/.grass7/addons/scripts/r.in.nasadem", line 342, in createTMPlocation if grass.parse_command('g.proj', flags='g')['EPSG'] != str(epsg): KeyError: 'EPSG' ERROR: Region <r_in_nasadem_region_409609> not found Error in atexit._run_exitfuncs: Traceback (most recent call last): File "/home/mundialis/.grass7/addons/scripts/r.in.nasadem", line 312, in cleanup grass.run_command("g.region", region=tmpregionname) File "/home/mundialis/software/grass78_git/dist.x86_64-pc-linux-gnu/etc/python/grass/script/core.py", line 441, in run_command return handle_errors(returncode, returncode, args, kwargs) File "/home/mundialis/software/grass78_git/dist.x86_64-pc-linux-gnu/etc/python/grass/script/core.py", line 342, in handle_errors raise CalledModuleError(module=None, code=code, grass.exceptions.CalledModuleError: Module run None g.region region=r_in_nasadem_region_409609 ended with error Process ended with non-zero return code 1. See errors in the (error) output. ``` The change is the introduction of `srid`: ``` {'name': 'WGS 84', 'datum': 'wgs84', 'ellps': 'wgs84', 'proj': 'll', 'no_defs': 'defined', 'srid': 'EPSG:4326', 'unit': 'degree', 'units': 'degrees', 'meters': '1.0'} ``` which is addressed by this PR (feedback concerning the Python style is welcome). Fixes - r.in.nasadem - r.in.srtm.region (perhaps other addons are affected as well) Co-authored-by: Guido Riembauer <62383722+griembauer@users.noreply.github.com>
Co-authored-by: nilason <n_larsson@yahoo.com> * adjust lib/gis, lib/proj, g.proj, r.proj, v.proj, r.in.gdal, r.external, r.out.gdal, v.in.ogr, v.external, v.out.ogr
Co-authored-by: nilason <n_larsson@yahoo.com> * adjust lib/gis, lib/proj, g.proj, r.proj, v.proj, r.in.gdal, r.external, r.out.gdal, v.in.ogr, v.external, v.out.ogr
New GRASS locations now store the projection definition also in WKT (WKT2 with GDAL 3 / PROJ 6) and the spatial reference ID (SRID) if available. Two new files are (if possible) created for a new location: PERMANENT/PROJ_SRID and PERMANENT/PROJ_WKT. The file PERMANENT/PROJ_EPSG is no longer created and superseeded by PROJ_SRID which supports not only EPSG, but also other authorities. An existing file PERMANENT/PROJ_EPSG is read if PROJ_SRID does not exist.
For reprojections, the priority of spatial reference definitions is 1. SRID, 2. WKT, 3. PROJ4 style from GRASS definition
changed libraries
gis, proj
changed modules
g.proj, r.proj, v.proj, r.in.gdal, v.in.ogr, r.external, v.external, r.out.gdal, v.out.ogr
This PR also solves qgis/QGIS#18596 (comment)
testing
reprojection from/to a CRS with neither meters nor degrees as units might not convert coordinates to map units
use location nc_spm_08
EPSG:3358
NAD83(HARN) / North Carolina
CRS with feet as units
EPSG:2264
NAD83 / North Carolina (ftUS)
location nc_epsg_2264
CRS with different datum
EPSG:32617
UTM zone 17N
location nc_epsg_32617
reprojection with projection definitions from SRIDs
reproject nc_state from nc_spm_08 to location nc_epsg_2264 to location nc_epsg_32617 and back to nc_spm_08
reproject nc_state from nc_spm_08 to location nc_epsg_32617 to location nc_epsg_2264 and back to nc_spm_08
end results in nc_spm_08 should be identical (within floating point precision limits)
now rename all PROJ_EPSG, PROJ_SRID, and PROJ_WKT files in the three locations
reprojection with standard GRASS projection definitions
reproject nc_state from nc_spm_08 to location nc_epsg_2264 to location nc_epsg_32617 and back to nc_spm_08
reproject nc_state from nc_spm_08 to location nc_epsg_32617 to location nc_epsg_2264 and back to nc_spm_08
end results in nc_spm_08 should be identical (within floating point precision limits)
This is not yet working properly because the PROJ area of use (proj_area_set_bbox(), proj_create_crs_to_crs()) is not yet set properly. The new PROJ API needs an area of use, similar to the GRASS current region, to select the best coordinate transformation method. If this area of use is not set correctly, a wrong coordinate transformation method might be selected, and reprojection can fail. Therefore this PR is WIP.