Context Navigation

← Previous Ticket
Next Ticket →

#3191 closed defect (wontfix)

coordinate reference system issue when exporting to Rasterlite

Reported by:	vincent	Owned by:	grass-dev@…
Priority:	normal	Milestone:	7.0.6
Component:	Raster	Version:	unspecified
Keywords:	gdal export, Rasterlite	Cc:
CPU:	Unspecified	Platform:	Unspecified

Description

r.out.gdal allows to export raster data to Rasterlite tables. It looks like the output raster table lacks a correct reference to the srid it should be linked to, i.e. in the geometry_column table, srid equals -1 for the newly written raster table.

Change History (8)

in reply to: description comment:1 by hellik, 8 years ago

Replying to vincent:

r.out.gdal allows to export raster data to Rasterlite tables. It looks like the output raster table lacks a correct reference to the srid it should be linked to, i.e. in the geometry_column table, srid equals -1 for the newly written raster table.

tested with

System Info                                                                     
GRASS version: 7.3.svn                                                          
GRASS SVN revision: r69714                                                      
Build date: 2016-10-24                                                          
Build platform: i386-w64-mingw32                                                
GDAL: 2.1.1                                                                     
PROJ.4: 4.9.3                                                                   
GEOS: 3.5.0                                                                     
SQLite: 3.14.1                                                                  
Python: 2.7.4                                                                   
wxPython: 2.8.12.1                                                              
Platform: Windows-Vista-6.0.6002-SP2 (OSGeo4W)

with data (downloadable here data.

-*--*--*--*--*--*--*--*--*--*-
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
-*--*--*--*--*--*--*--*--*--*-
projection: 99 (MGI / Austria GK Central)
zone:       0
datum:      hermannskogel
ellipsoid:  bessel
north:      230005
south:      159995
west:       -96005
east:       -23995
nsres:      10
ewres:      10
rows:       7001
cols:       7201
cells:      50414201
-*--*--*--*--*--*--*--*--*--*-
r.out.gdal input=lz_10m_float@user output=J:\wd\checkraster\lz_10m\lz_10m\geotiff_out.tif format=GTiff
-*--*--*--*--*--*--*--*--*--*-
gdalinfo geotiff_out.tif
Driver: GTiff/GeoTIFF
Files: geotiff_out.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)
Band 1 Block=7201x1 Type=Float32, ColorInterp=Gray
  NoData Value=nan
-*--*--*--*--*--*--*--*--*--*-
gdal_translate -of Rasterlite geotiff_out.tif gdtrans_to_rasterlite.db
-*--*--*--*--*--*--*--*--*--*-
gdalinfo gdtrans_to_rasterlite.db
Driver: Rasterlite/Rasterlite
Files: gdtrans_to_rasterlite.db
Size is 7201, 7001
Coordinate System is `'
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  TILE_FORMAT=GTiff
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -96005.000,  230005.000)
Lower Left  (  -96005.000,  159995.000)
Upper Right (  -23995.000,  230005.000)
Lower Right (  -23995.000,  159995.000)
Center      (  -60000.000,  195000.000)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
-*--*--*--*--*--*--*--*--*--*-
gdal_translate -a_srs EPSG:31255 -of Rasterlite geotiff_out.tif gdtrans_to_rasterlite_srs.db
-*--*--*--*--*--*--*--*--*--*-
gdalinfo gdtrans_to_rasterlite_srs.db
Driver: Rasterlite/Rasterlite
Files: gdtrans_to_rasterlite_srs.db
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"]],
        AXIS["Latitude",NORTH],
        AXIS["Longitude",EAST],
        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"]],
    AXIS["X",NORTH],
    AXIS["Y",EAST],
    AUTHORITY["EPSG","31255"]]
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  TILE_FORMAT=GTiff
Image Structure Metadata:
  INTERLEAVE=PIXEL
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)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
-*--*--*--*--*--*--*--*--*--*-
r.out.gdal input=lz_10m_float@user output=grass_to_rasterlite2.db format=Rasterlite
-*--*--*--*--*--*--*--*--*--*-
gdalinfo grass_to_rasterlite2.db
Driver: Rasterlite/Rasterlite
Files: grass_to_rasterlite2.db
Size is 7201, 7001
Coordinate System is `'
Origin = (-96005.000000000000000,230005.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  TILE_FORMAT=GTiff
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -96005.000,  230005.000)
Lower Left  (  -96005.000,  159995.000)
Upper Right (  -23995.000,  230005.000)
Lower Right (  -23995.000,  159995.000)
Center      (  -60000.000,  195000.000)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
-*--*--*--*--*--*--*--*--*--*-

testing with gdal_translate and r.out.gdal shows

a simple gdal_translate from geotiff to rasterlite without assigning the srs, the resulting rasterlite db has no srs info

gdal_translate with e.g. -a_srs EPSG:your_epsg_code, the resulting rasterlite db has srs info.

so it seems to be a gdal "issue"; maybe it's also worth to ask in the GDAL ML?

follow-up: 4 comment:2 by vincent, 8 years ago

Yes, I've come to the same observation so far, but I don't know enough about how GRASS parses projection info to gdal; sorry for this rather naive question but in the case of GTiff driver why then is there no need to worry about coordinate reference system information for it to work?

Well, if you all prefer, we can close the ticket.

comment:3 by sbl, 8 years ago

Maybe related to: #3059
Symptoms are the same...

in reply to: 2 ; follow-up: 6 comment:4 by hellik, 8 years ago

Replying to vincent:

Yes, I've come to the same observation so far, but I don't know enough about how GRASS parses projection info to gdal; sorry for this rather naive question but in the case of GTiff driver why then is there no need to worry about coordinate reference system information for it to work?

Well, if you all prefer, we can close the ticket.

taken from the GDAL ML

> is this intentional that there is no srs by a simple gdal_translate of a
> geotiff to a Rasterlite db?

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. 
From what I can see, it also probably doesn't properly support creating 
entries in spatial_ref_sys with Spatialite 4 layout (the driver has been 
likely developped during Spatialite 3 era)
It is "intentional" in the meaning that this wasn't implemented, but could 
potentially be added. Could be worth a ticket in Trac.

> maybe there is some other command line option related to Rasterlite db
> which I'm missing?

No, the workaround is the one you found. Make sure a EPSG code is attached to 
the input dataset, or assign it explictly

I think this ticket can be closed; but I'll open another one for the last issue mentioned above.

comment:5 by hellik, 8 years ago

Resolution:	→ wontfix
Status:	new → closed

in reply to: 4 ; follow-up: 7 comment:6 by hellik, 8 years ago

Replying to hellik:

I think this ticket can be closed; but I'll open another one for the last issue mentioned above.

ticket opened; see #3193

in reply to: 6 ; follow-up: 8 comment:7 by vincent, 8 years ago

Replying to hellik:

Replying to hellik:

I think this ticket can be closed; but I'll open another one for the last issue mentioned above.

ticket opened; see #3193

Not sure I understand: do you mean adding this "AUTHORITY node" ability to r.out.gdal may help in the present case? or is there really something wrong with gdal?

in reply to: 7 comment:8 by hellik, 8 years ago

Replying to vincent:

Replying to hellik:

Replying to hellik:

I think this ticket can be closed; but I'll open another one for the last issue mentioned above.

ticket opened; see #3193

Not sure I understand: do you mean adding this "AUTHORITY node" ability to r.out.gdal may help in the present case? or is there really something wrong with gdal?

regarding "AUTHORITY node", maybe yes.

as Even mentioned, there is nothing wrong in gdal, just not implemented.

Note: See TracTickets for help on using tickets.

Download in other formats: