Opened 3 years ago
Closed 15 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 23 months ago by
This functionality was recently implemented in trunk, right? Reference is missing here.
comment:6 follow-up: 9 Changed 23 months ago by
comment:7 Changed 23 months 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 23 months ago by
This affects interoperability with e.g. GeoServer/GeoNode? or QGIS (https://issues.qgis.org/issues/14856).
comment:9 follow-up: 10 Changed 23 months 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 23 months 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:11 Changed 23 months ago by
Milestone: | 7.4.0 → 7.4.1 |
---|
Ticket retargeted after milestone closed
comment:12 follow-up: 13 Changed 23 months 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 23 months 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 22 months 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 18 months ago by
Milestone: | 7.4.1 → 7.4.2 |
---|
comment:16 Changed 15 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