[GRASS-user] Grass7 - Inconsistencies with v.out.ogr and v.out.postgis

Possibly some bugs. In summary:

CASE 1: v.out.ogr if I export to shapefile works as expected - in terms of 176783 features and SRID 2193.
CASE 2: v.out.ogr if I export to postgis, only returns 139 features, with SRID 900915(?)
CASE 3: v.out.postgis returns 176782 features, with SRID=0 despite setting SRID=2193 using the options flag.

Has anyone else encountered similar problems?

The detail...

-----------
CASE 1
-----------
v.category osm_nz_roads_2193_cleaned op=report
Layer/table: 1/grass.osm_nz_roads_2193_cleaned
type count min max
point 0 0 0
line 176783 1 88617
boundary 0 0 0
centroid 0 0 0
area 0 0 0
face 0 0 0
kernel 0 0 0
all 176783 1 88617

v.out.ogr --overwrite input=osm_nz_roads_2193_cleaned type=line dsn="PG:host=localhost dbname=test user=postgres password=password active_schema=osm" olayer=osm_nz_roads_2193_export format=PostgreSQL

Rows (estimated) 139

-- Table: osm.osm_nz_roads_2193_export

-- DROP TABLE osm.osm_nz_roads_2193_export;

CREATE TABLE osm.osm_nz_roads_2193_export
(
  ogc_fid serial NOT NULL,
  wkb_geometry geometry(LineString,900915),
  cat integer,
  osm_id integer,
  name character varying(255),
  highway character varying(255),
  bridge integer,
  oneway integer,
  driveable integer,
  trainable integer,
  walkable integer,
  cycleable integer,
  ignore integer,
  CONSTRAINT osm_nz_roads_2193_export_pk PRIMARY KEY (ogc_fid)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE osm.osm_nz_roads_2193_export
  OWNER TO postgres;

-- Index: osm.osm_nz_roads_2193_export_geom_idx

-- DROP INDEX osm.osm_nz_roads_2193_export_geom_idx;

CREATE INDEX osm_nz_roads_2193_export_geom_idx
  ON osm.osm_nz_roads_2193_export
  USING gist
  (wkb_geometry);

-----------
CASE 2
-----------

v.out.postgis --o input=osm_nz_roads_2193_cleaned dsn="PG:host=localhost dbname=test user=postgres password=password" olayer=osm.osm_nz_roads_2193_export options="GEOMETRY_NAME=wkb_geometry,SRID=2193"

Rows (estimated) 176782

- Table: osm.osm_nz_roads_2193_export

-- DROP TABLE osm.osm_nz_roads_2193_export;

CREATE TABLE osm.osm_nz_roads_2193_export
(
  fid serial NOT NULL,
  cat integer,
  osm_id integer,
  name character varying(255),
  highway character varying(255),
  bridge integer,
  oneway integer,
  driveable integer,
  trainable integer,
  walkable integer,
  cycleable integer,
  ignore integer,
  wkb_geometry geometry(LineString),
  CONSTRAINT osm_nz_roads_2193_export_pkey PRIMARY KEY (fid)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE osm.osm_nz_roads_2193_export
  OWNER TO postgres;

-- Index: osm.osm_nz_roads_2193_export_wkb_geometry_idx

-- DROP INDEX osm.osm_nz_roads_2193_export_wkb_geometry_idx;

CREATE INDEX osm_nz_roads_2193_export_wkb_geometry_idx
  ON osm.osm_nz_roads_2193_export
  USING gist
  (wkb_geometry);

-----------
CASE 3
-----------
v.out.ogr input=osm_nz_roads_2193_cleaned type=line dsn=/var/tmp/osm_nz_roads_export.shp

ogrinfo -so osm_nz_roads_export.shp osm_nz_roads_export
INFO: Open of `osm_nz_roads_export.shp'
      using driver `ESRI Shapefile' successful.

Layer name: osm_nz_roads_export
Geometry: Line String
Feature Count: 176783
Extent: (1148817.072845, 4793789.625743) - (2089039.165212, 6193819.546965)
Layer SRS WKT:
PROJCS["Transverse_Mercator",
    GEOGCS["GCS_grs80",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["Meter",1]]
cat: Integer (10.0)
osm_id: Integer (10.0)
name: String (255.0)
highway: String (255.0)
bridge: Integer (10.0)
oneway: Integer (10.0)
driveable: Integer (10.0)
trainable: Integer (10.0)
walkable: Integer (10.0)
cycleable: Integer (10.0)
ignore: Integer (10.0)

Hi,

2013/4/26 Mark Wynter <mark@dimensionaledge.com>:

Possibly some bugs. In summary:

CASE 1: v.out.ogr if I export to shapefile works as expected - in terms of 176783 features and SRID 2193.
CASE 2: v.out.ogr if I export to postgis, only returns 139 features, with SRID 900915(?)
CASE 3: v.out.postgis returns 176782 features, with SRID=0 despite setting SRID=2193 using the options flag.

Has anyone else encountered similar problems?

note that `v.out.postgis` is quite new module (mostly related to
ongoing PostGIS Topology support). Generally speaking work in
progress. I will look at these issues ASAP. Thanks for testing. Martin

Hi,

2013/4/26 Mark Wynter <mark@dimensionaledge.com>:

v.category osm_nz_roads_2193_cleaned op=report
Layer/table: 1/grass.osm_nz_roads_2193_cleaned
type count min max
point 0 0 0
line 176783 1 88617
boundary 0 0 0
centroid 0 0 0
area 0 0 0
face 0 0 0
kernel 0 0 0
all 176783 1 88617

note that `v.out.ogr` export by default only features with category
(see `-c` flag).

I have tried `roadsmajor` from `nc` dataset.

$ v.category roadsmajor opt=r -g
1 line 355 1 355
1 all 355 1 355

v.out.ogr --overwrite input=osm_nz_roads_2193_cleaned type=line dsn="PG:host=localhost dbname=test user=postgres password=password active_schema=osm" olayer=osm_nz_roads_2193_export format=PostgreSQL

v.out.ogr in=roadsmajor dsn=PG:dbname=pgis_grass format=PostgreSQL

->

$ db.select data=pgis_grass dri=pg sql="select count(*) from roadsmajor" -c
355

db.select data=pgis_grass dri=pg sql="select srid from
geometry_columns where f_table_name = 'roadsmajor'" -c
900914

(new srid added to the table `spatial_ref_sys`)

v.out.postgis --o input=osm_nz_roads_2193_cleaned dsn="PG:host=localhost dbname=test user=postgres password=password" olayer=osm.osm_nz_roads_2193_export options="GEOMETRY_NAME=wkb_geometry,SRID=2193"

$ v.out.postgis in=roadsmajor dsn=PG:dbname=pgis_grass
options="GEOMETRY_NAME=wkb_geometry,SRID=3358" --o

$ db.select data=pgis_grass dri=pg sql="select count(*) from roadsmajor" -c
355

$ db.select data=pgis_grass dri=pg sql="select srid from
geometry_columns where f_table_name = 'roadsmajor'" -c
3358

v.out.ogr input=osm_nz_roads_2193_cleaned type=line dsn=/var/tmp/osm_nz_roads_export.shp

$ v.out.ogr in=roadsmajor dsn=/tmp/roadsmajor.shp

ogrinfo -so osm_nz_roads_export.shp osm_nz_roads_export

$ ogrinfo /tmp/roadsmajor.shp -so roadsmajor
INFO: Open of `/tmp/roadsmajor.shp'
      using driver `ESRI Shapefile' successful.

Layer name: roadsmajor
Geometry: Line String
Feature Count: 355
Extent: (611136.499873, 197465.257810) - (676800.487681, 257970.129540)
Layer SRS WKT:
PROJCS["Lambert_Conformal_Conic",
    GEOGCS["GCS_grs80",
        DATUM["North_American_Datum_1983",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",36.16666666666666],
    PARAMETER["standard_parallel_2",34.33333333333334],
    PARAMETER["latitude_of_origin",33.75],
    PARAMETER["central_meridian",-79],
    PARAMETER["false_easting",609601.22],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
cat: Integer (10.0)
MAJORRDS_: Real (24.15)
ROAD_NAME: String (16.0)
MULTILANE: String (3.0)
PROPYEAR: Integer (10.0)
OBJECTID: Integer (10.0)
SHAPE_LEN: Real (24.15)

Everything seems to be OK. Do you have fresh installation from svn? Martin

--
Martin Landa <landa.martin gmail.com> * http://geo.fsv.cvut.cz/~landa

Do you have fresh installation from svn?

Thanks Martin. Good advice - I just pulled the latest from svn!

v.out.postgis now works, generating correct feature count and SRID=2193

yet... v.out.ogr when writing directly to PostGIS assigns SRID 900915 to the geometry column - when it should be 2193.
SRID 2193 is in the PostGIS spatial_ref_sys table.

v.out.ogr when writing to shapefile correctly assigns 2193 - see below.

Any idea as to why this is happening with v.out.ogr and PostGIS specifically? I suspect it has something to do with SRID 2193 not being recognised (somewhere in the process), hence a new SRID is created? Yet it works with v.out.postgis.

SELECT postgis_full_version();

                                                                           postgis_full_version
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
POSTGIS="2.1.0SVN r11004" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER
(1 row)

v.out.ogr input=osm_nz_roads_2193_cleaned type=line dsn=/var/tmp/osm_nz_roads_export.shp

ogrinfo -so osm_nz_roads_export.shp osm_nz_roads_export
INFO: Open of `osm_nz_roads_export.shp'
     using driver `ESRI Shapefile' successful.

Layer name: osm_nz_roads_export
Geometry: Line String
Feature Count: 176783
Extent: (1148817.072845, 4793789.625743) - (2089039.165212, 6193819.546965)
Layer SRS WKT:
PROJCS["Transverse_Mercator",
   GEOGCS["GCS_grs80",
       DATUM["New_Zealand_Geodetic_Datum_2000",
           SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
       PRIMEM["Greenwich",0],
       UNIT["Degree",0.017453292519943295]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["latitude_of_origin",0],
   PARAMETER["central_meridian",173],
   PARAMETER["scale_factor",0.9996],
   PARAMETER["false_easting",1600000],
   PARAMETER["false_northing",10000000],
   UNIT["Meter",1]]
cat: Integer (10.0)
osm_id: Integer (10.0)
name: String (255.0)
highway: String (255.0)
bridge: Integer (10.0)
oneway: Integer (10.0)
driveable: Integer (10.0)
trainable: Integer (10.0)
walkable: Integer (10.0)
cycleable: Integer (10.0)
ignore: Integer (10.0)

v.out.ogr in=roadsmajor dsn=PG:dbname=pgis_grass format=PostgreSQL

$ db.select data=pgis_grass dri=pg sql=“select count(*) from roadsmajor” -c
355

db.select data=pgis_grass dri=pg sql=“select srid from
geometry_columns where f_table_name = ‘roadsmajor’” -c
900914

(new srid added to the table spatial_ref_sys)

Hi Martin, just noticed PostGIS assigned 900914 to your table above, yet below its 3358.
This is the same as what I’m experiencing.
When I inspect the details of 900914, it mirrors the correct SRID - its just labelled 900914 instead?
Yet v.out.postgis handles everything the way you would expect

$ v.out.postgis in=roadsmajor dsn=PG:dbname=pgis_grass
options=“GEOMETRY_NAME=wkb_geometry,SRID=3358” --o

$ db.select data=pgis_grass dri=pg sql=“select count(*) from roadsmajor” -c
355

$ db.select data=pgis_grass dri=pg sql=“select srid from
geometry_columns where f_table_name = ‘roadsmajor’” -c
3358

Everything seems to be OK. Do you have fresh installation from svn? Martin

Mark wrote:

-- Table: osm.osm_nz_roads_2193_export

(putting on another hat)

by the way, for NZ roads in OSM we have the full set of updated
gov't road data (cleaned and connected by the NZOpenGPS group)
sitting ready for merge. So far we've been concentrating on
uploading the natural layers to OSM as they require less merging
work with existing data, but all are welcome to work on whatever
they please.

see http://linz2osm.openstreetmap.org.nz/ for the latest data
and status, grouped by provincial region.

There is a several years old version of the NZ road data here:
  http://joerichards.dev.openstreetmap.org/files.html

but it is complete and ready for loading into PostGIS with
osm2pgsql.

more info at http://wiki.openstreetmap.org/wiki/LINZ
coordination of efforts is done on the NZOpenGIS google group.

regards,
Hamish

Hi,

2013/4/27 Mark Wynter <mark@dimensionaledge.com>:

Do you have fresh installation from svn?

Thanks Martin. Good advice - I just pulled the latest from svn!

v.out.postgis now works, generating correct feature count and SRID=2193

the source of problem is missing EPSG in your location. There was
discussion about that in ML AFAIR some months ago, but with no result.
`v.out.ogr` knows which projection you are using, but doesn't know
which EPSG it could be. So when exporting data to PostGIS using
`v.out.ogr` the tries to find matching projection in
`spatial_ref_sys`. If no matching projection is found, then the new
projection is added to `spatial_ref_sys`.

Eg. when working in `nc` dataset (epsg: 3358)

$ db.select data=pgis_grass dri=pg sql="select proj4text from
spatial_ref_sys where srid = 3358" -c
+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334
+lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs

$ v.out.ogr in=roadsmajor dsn=PG:dbname=pgis_grass format=PostgreSQL

Current projection is:

g.proj -jf
+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334
+lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +no_defs +a=6378137
+rf=298.257222101 +towgs84=0.000,0.000,0.000 +to_meter=1

There is no exact match between them, so new srid is added to `spatial_ref_sys`.

$ db.select data=pgis_grass dri=pg sql="select proj4text from
spatial_ref_sys where srid = 900914" -c
+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334
+lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m
+no_defs

yet... v.out.ogr when writing directly to PostGIS assigns SRID 900915 to the geometry column - when it should be 2193.
SRID 2193 is in the PostGIS spatial_ref_sys table.

`v.out.postgis` behaves in similar way, unless you define `options=srid=3358`.

v.out.ogr when writing to shapefile correctly assigns 2193 - see below.

? Naturally there is no information about EPSG.

ogrinfo /tmp/roadsmajor.shp -so roadsmajor
INFO: Open of `/tmp/roadsmajor.shp'
      using driver `ESRI Shapefile' successful.

Layer name: roadsmajor
Geometry: Line String
Feature Count: 355
Extent: (611136.499873, 197465.257810) - (676800.487681, 257970.129540)
Layer SRS WKT:
PROJCS["Lambert_Conformal_Conic",
    GEOGCS["GCS_grs80",
        DATUM["North_American_Datum_1983",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",36.16666666666666],
    PARAMETER["standard_parallel_2",34.33333333333334],
    PARAMETER["latitude_of_origin",33.75],
    PARAMETER["central_meridian",-79],
    PARAMETER["false_easting",609601.22],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
cat: Integer (10.0)
MAJORRDS_: Real (24.15)
ROAD_NAME: String (16.0)
MULTILANE: String (3.0)
PROPYEAR: Integer (10.0)
OBJECTID: Integer (10.0)
SHAPE_LEN: Real (24.15)

Martin

--
Martin Landa <landa.martin gmail.com> * http://geo.fsv.cvut.cz/~landa

Martin wrote:

the source of problem is missing EPSG in your location.

I'd suggest the problem is rather in the programming assumptions.

There was discussion about that in ML AFAIR some months ago,
but with no result.

sure, the result was if a location is created with an EPSG code, the
EPSG code is recorded into a PERMANENT/PROJ_EPSG file. That file is
not, and should not, be used for anything other than historical purposes. (it's a one way street, grass defines the projection
terms with more precision than the EPSG codes does [e.g. transform terms
used])

That file should only be used for historical metadata
/user interest pruposes. All projection code should go through
the pj_*() functions.

`v.out.ogr` knows which projection you are using, but doesn't know
which EPSG it could be. So when exporting data to PostGIS using
`v.out.ogr` the tries to find matching projection in
`spatial_ref_sys`.

that is a flawed approach, it's a near impossible task both in process
and in result. See the fun that happened some years ago when Frank
changed the precision of the proj4 epsg file by a significant digit
or two, everyone who did string matching broke, those which tested
abs(dx)<epsilon or some sort (if there were any) found out if their
epsilon was set right etc.
You'd have to know about the full datum -> ellipsoid+transform
terms, and back (but you can't go back because the transform terms
are local and not formally defined). And even then, order of the terms
matter, so you can't just sort them and test term by term. While
not impossible, it is just such an ugly road to go down the best
thing is to not try rather than set ourselves up to fail and end
up with almost certainly bug riddled software.

If no matching projection is found, then the new projection is
added to `spatial_ref_sys`.

eek! that's an even worse idea. Those are definitive codes, nobody but
the EPSG people should be assigning them. We are dealing with formal
definitions and we are not a granting body. AFAIK PostGIS can take
a full proj4 string to define the SRS, so why start making up EPSG
codes? Just figure out the way to use the full string. (which
isn't perfect either, as GRASS's projection def'ns can contain
more detail than +proj terms know how to express)
Defining new SRS codes as anything more than a quick hack is
simply wrong.

trying to save some later grief,
Hamish