[GRASS-dev] [GRASS GIS] #3191: coordinate reference system issue when exporting to Rasterlite

GRASS GIS trac at osgeo.org
Mon Oct 24 04:54:52 PDT 2016


#3191: coordinate reference system issue when exporting to Rasterlite
--------------------------+-------------------------------------
  Reporter:  vincent      |      Owner:  grass-dev@…
      Type:  defect       |     Status:  new
  Priority:  normal       |  Milestone:  7.0.6
 Component:  Raster       |    Version:  unspecified
Resolution:               |   Keywords:  gdal export, Rasterlite
       CPU:  Unspecified  |   Platform:  Unspecified
--------------------------+-------------------------------------

Comment (by hellik):

 Replying to [ticket:3191 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 [https://www.tirol.gv.at/data/datenkatalog
 /geographie-und-planung/digitales-gelaendemodell-tirol/ 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 at 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 at 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?

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3191#comment:1>
GRASS GIS <https://grass.osgeo.org>



More information about the grass-dev mailing list