Opened 3 years ago
Closed 16 months ago
#3193 closed enhancement (fixed)
r.out.gdal: add AUTHORITY node to the srs info of the exported raster
| Reported by: | hellik | Owned by: | |
|---|---|---|---|
| Priority: | major | Milestone: | 7.6.0 |
| Component: | Raster | Version: | svn-releasebranch76 |
| Keywords: | r.external.out, r.out.gdal, v.external.out, v.out.ogr | Cc: | |
| CPU: | All | Platform: | All |
Description
test case in a location created by an EPSG code
g.proj -p -PROJ_INFO------------------------------------------------- name : MGI / Austria GK Central datum : hermannskogel ellps : bessel proj : tmerc lat_0 : 0 lon_0 : 13.33333333333333 k : 1 x_0 : 0 y_0 : -5000000 no_defs : defined towgs84 : 577.326,90.129,463.919,5.1366,1.4742,5.2970,2.4232 -PROJ_EPSG------------------------------------------------- epsg : 31255 -PROJ_UNITS------------------------------------------------ unit : meter units : meters meters : 1
export a raster by r.out.gdal and do gdalinfo
gdalinfo raster_export_by_grass.tif
Driver: GTiff/GeoTIFF
Files: raster_export_by_grass.tif
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
GEOGCS["bessel",
DATUM["Militar_Geographische_Institute",
SPHEROID["Bessel_1841",6377397.155,299.1528128],
TOWGS84[577.326,90.129,463.919,5.1366,1.4742,5.297,2.4232]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",13.33333333333333],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",-5000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=GRASS GIS 7.3.svn with GDAL 2.1.1
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -96005.000, 230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left ( -96005.000, 159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right ( -23995.000, 230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right ( -23995.000, 159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center ( -60000.000, 195000.000) ( 12d32'45.35"E, 46d53'30.77"N)
do gdal_translate and assign the srs
gdal_translate -a_srs EPSG:31255 raster_export_by_grass.tif geotif_by_gdaltranslate.tif
check by gdalinfo
gdalinfo geotif_by_gdaltranslate.tif
Driver: GTiff/GeoTIFF
Files: geotif_by_gdaltranslate.tif
geotif_by_gdaltranslate.tif.aux.xml
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
GEOGCS["MGI",
DATUM["Militar_Geographische_Institute",
SPHEROID["Bessel 1841",6377397.155,299.1528128,
AUTHORITY["EPSG","7004"]],
TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],
AUTHORITY["EPSG","6312"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4312"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",13.33333333333333],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",-5000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","31255"]] <<==
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=GRASS GIS 7.3.svn with GDAL 2.1.1
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -96005.000, 230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left ( -96005.000, 159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right ( -23995.000, 230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right ( -23995.000, 159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center ( -60000.000, 195000.000) ( 12d32'45.35"E, 46d53'30.77"N)
GDAL adds the AUTHORITY node of the corresponding EPSG code.
this may be help #3191 (and maybe also other GIS software reading GRASS exported data).
also mentioned in a GDAL ML thread:
From what I can see the driver doesn't handle adding a new entry in the spatial_ref_sys table if there's no AUTHORITY node, such as in your use case. [...] Make sure a EPSG code is attached to the input dataset, or assign it explictly
Change History (16)
comment:1 Changed 3 years ago by
| Keywords: | r.out.gdal r.external.out added |
|---|
comment:3 Changed 3 years ago by
| Milestone: | 7.2.1 → 7.2.2 |
|---|
comment:4 Changed 2 years ago by
| Milestone: | 7.2.2 → 7.4.0 |
|---|
All enhancement tickets should be assigned to 7.4 milestone.
comment:5 follow-up: 6 Changed 2 years ago by
This functionality was recently implemented in trunk, right? Reference is missing here.
comment:6 follow-up: 9 Changed 2 years ago by
comment:7 Changed 2 years ago by
Hmm... This is still an issue for me both in GRASS 7.4 (r72108) and GRASS 7.5 (r72157):
See also (sort of duplicate): #3048
g.proj -p -PROJ_INFO------------------------------------------------- name : ETRS89 / UTM zone 33N datum : etrs89 ellps : grs80 proj : utm zone : 33 no_defs : defined -PROJ_EPSG------------------------------------------------- epsg : 25833 -PROJ_UNITS------------------------------------------------ unit : meter units : meters meters : 1 r.external.out --verbose directory=/tmp extension=.tif format=GTiff options=COMPRESS=LZW r.mapcalc expression="test=1" gdalinfo /tmp/test.tif r.external.out -r g.remove type=raster name=test -f rm /tmp/test.tif r.mapcalc expression="test=1" r.out.gdal -c input=test output=/tmp/test.tif createopt="COMPRESS=LZW" --o gdalinfo /tmp/test.tif
In both cases gdalinfo shows:
Coordinate System is:
PROJCS["UTM Zone 33, Northern Hemisphere",
GEOGCS["grs80",
DATUM["European_Terrestrial_Reference_System_1989",
SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6258"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",15],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (256495.000000000000000,6653695.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
That means:
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","25833"]]
is still missing at the bottom of the CRS information, compared to CRS info with:
gdal_edit.py -a_srs EPSG:25833 /tmp/test.tif gdalinfo /tmp/test.tif
comment:8 Changed 2 years ago by
This affects interoperability with e.g. GeoServer/GeoNode? or QGIS (https://issues.qgis.org/issues/14856).
comment:9 follow-up: 10 Changed 2 years ago by
Replying to mmetz:
Replying to martinl:
This functionality was recently implemented in trunk, right? Reference is missing here.
It was implemented by martinl with r71303
tested with the example above:
g.proj -p -PROJ_INFO------------------------------------------------- name : MGI / Austria GK Central datum : hermannskogel ellps : bessel proj : tmerc lat_0 : 0 lon_0 : 13.33333333333333 k : 1 x_0 : 0 y_0 : -5000000 no_defs : defined towgs84 : 577.326,90.129,463.919,5.1366,1.4742,5.2970,2.4232 -PROJ_EPSG------------------------------------------------- epsg : 31255 -PROJ_UNITS------------------------------------------------ unit : meter units : meters meters : 1
r.out.gdal input=test1@test output=J:\test\test_write_auth.tif format=GTiff createopt=COMPRESS=LZW,PREDICTOR=2,TILED=YES
gdalinfo test_write_auth.tif
Driver: GTiff/GeoTIFF
Files: test_write_auth.tif
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
GEOGCS["bessel",
DATUM["Militar_Geographische_Institute",
SPHEROID["Bessel_1841",6377397.155,299.1528128],
TOWGS84[577.326,90.129,463.919,5.1366,1.4742,5.297,2.4232]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",13.33333333333333],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",-5000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=GRASS GIS 7.5.svn with GDAL 2.2.3
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -96005.000, 230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left ( -96005.000, 159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right ( -23995.000, 230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right ( -23995.000, 159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center ( -60000.000, 195000.000) ( 12d32'45.35"E, 46d53'30.77"N)
gdal_translate -a_srs EPSG:31255 test_write_auth.tif test_write_auth_assign_by_gdal.tif
Driver: GTiff/GeoTIFF
Files: test_write_auth_assign_by_gdal.tif
Size is 7201, 7001
Coordinate System is:
PROJCS["MGI / Austria GK Central",
GEOGCS["MGI",
DATUM["Militar_Geographische_Institute",
SPHEROID["Bessel 1841",6377397.155,299.1528128,
AUTHORITY["EPSG","7004"]],
TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],
AUTHORITY["EPSG","6312"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4312"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",13.33333333333333],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",-5000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","31255"]] <==
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=GRASS GIS 7.5.svn with GDAL 2.2.3
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -96005.000, 230005.000) ( 12d 3'57.90"E, 47d12' 8.96"N)
Lower Left ( -96005.000, 159995.000) ( 12d 4'50.98"E, 46d34'22.11"N)
Upper Right ( -23995.000, 230005.000) ( 13d 0'59.64"E, 47d12'32.62"N)
Lower Right ( -23995.000, 159995.000) ( 13d 1'12.91"E, 46d34'45.25"N)
Center ( -60000.000, 195000.000) ( 12d32'45.35"E, 46d53'30.77"N)
System Info GRASS version: 7.5.svn GRASS SVN revision: r72120 Build date: 2018-01-24 Build platform: i386-w64-mingw32 GDAL: 2.2.3 PROJ.4: 4.9.3 GEOS: 3.5.0 SQLite: 3.17.0 Python: 2.7.4 wxPython: 2.8.12.1 Platform: Windows-Vista-6.0.6002-SP2 (OSGeo4W)
comparing the gdalinfo of GRASS exported and GDAL assigned tif,
AUTHORITY["EPSG","31255"]
seems to be still missing
comment:10 follow-ups: 12 14 Changed 2 years ago by
Replying to hellik:
Replying to mmetz:
Replying to martinl:
This functionality was recently implemented in trunk, right? Reference is missing here.
It was implemented by martinl with r71303
tested with the example above: [...]
AUTHORITY["EPSG","31255"]seems to be still missing
The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.
comment:12 follow-up: 13 Changed 2 years ago by
Replying to mmetz:
Replying to hellik:
Replying to mmetz:
Replying to martinl:
This functionality was recently implemented in trunk, right? Reference is missing here.
It was implemented by martinl with r71303
tested with the example above: [...]
AUTHORITY["EPSG","31255"]seems to be still missing
The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.
Backport candidat for 7.4.1?
comment:13 Changed 2 years ago by
Replying to hellik:
Replying to mmetz:
Replying to hellik:
Replying to mmetz:
Replying to martinl:
This functionality was recently implemented in trunk, right? Reference is missing here.
It was implemented by martinl with r71303
tested with the example above: [...]
AUTHORITY["EPSG","31255"]seems to be still missing
The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.
Backport candidat for 7.4.1?
It's a new feature, and EPSG codes should not be required in order to support custom coordinate reference systems.
comment:14 Changed 2 years ago by
Replying to mmetz:
Replying to hellik:
Replying to mmetz:
Replying to martinl:
This functionality was recently implemented in trunk, right? Reference is missing here.
It was implemented by martinl with r71303
tested with the example above: [...]
AUTHORITY["EPSG","31255"]seems to be still missing
The functions were implemented in r71303. An EPSG code was only exported with v.out.ogr. As of trunk r72160, an existing EPSG code is exported with r.external.out, r.out.gdal, v.external.out, v.out.ogr.
tested with
System Info GRASS version: 7.5.svn GRASS SVN revision: r72194 Build date: 2018-01-30 Build platform: x86_64-w64-mingw32 GDAL: 2.2.3 PROJ.4: 4.9.3 GEOS: 3.5.0 SQLite: 3.17.0 Python: 2.7.5 wxPython: 2.8.12.1 Platform: Windows-8-6.2.9200 (OSGeo4W)
g.proj -p -PROJ_INFO------------------------------------------------- name : MGI / Austria GK Central datum : hermannskogel ellps : bessel proj : tmerc lat_0 : 0 lon_0 : 13.33333333333333 k : 1 x_0 : 0 y_0 : -5000000 no_defs : defined towgs84 : 577.326,90.129,463.919,5.1366,1.4742,5.2970,2.4232 -PROJ_EPSG------------------------------------------------- epsg : 31255 -PROJ_UNITS------------------------------------------------ unit : meter units : meters meters : 1
gdalinfo test_raster_authority.tif
Driver: GTiff/GeoTIFF
Files: test_raster_authority.tif
Size is 170, 141
Coordinate System is:
PROJCS["MGI / Austria GK Central",
GEOGCS["MGI",
DATUM["Militar_Geographische_Institute",
SPHEROID["Bessel 1841",6377397.155,299.1528128,
AUTHORITY["EPSG","7004"]],
TOWGS84[577.326,90.129,463.919,5.1366,1.4742,5.297,2.4232],
AUTHORITY["EPSG","6312"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4312"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",13.33333333333333],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",-5000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","31255"]] <=
Origin = (-161500.000000000000000,273250.000000000000000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=GRASS GIS 7.5.svn with GDAL 2.2.3
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -161500.000, 273250.000) ( 11d11'10.74"E, 47d34'42.28"N)
Lower Left ( -161500.000, 238000.000) ( 11d11'56.99"E, 47d15'41.55"N)
Upper Right ( -119000.000, 273250.000) ( 11d45' 3.81"E, 47d35'15.35"N)
Lower Right ( -119000.000, 238000.000) ( 11d45'37.92"E, 47d16'14.25"N)
Center ( -140250.000, 255625.000) ( 11d28'27.33"E, 47d25'29.61"N)
AUTHORITY["EPSG","XXXX"] is now exported in rasters too.
comment:15 Changed 20 months ago by
| Milestone: | 7.4.1 → 7.4.2 |
|---|
comment:16 Changed 16 months ago by
| Keywords: | v.external.out v.out.ogr added |
|---|---|
| Milestone: | 7.4.2 → 7.6.0 |
| Resolution: | → fixed |
| Status: | new → closed |
| Version: | svn-trunk → svn-releasebranch76 |
To be released in the upcoming 7.6.0; closing.
See also: #3059