Opened 7 years ago

Closed 6 years ago

#3414 closed defect (fixed)

v.import: ERROR: Invalid e-w resolution field: 0 when importing EU_LAEA point into Sinusoidal location

Reported by: neteler Owned by: grass-dev@…
Priority: normal Milestone: 7.4.0
Component: Projections/Datums Version: svn-releasebranch74
Keywords: v.import Cc:
CPU: Unspecified Platform: Linux

Description

When trying to import a vector point in EU LAEA projection (see attached SHP file) into a location in MODIS Sinusoidal projection, an error occurs while v.proj does the job correctly:

GRASS 7.3.svn (modis_sinusoidal):~ > v.import in=point_eu_laea.shp out=point_eu_laea2
ERROR: Invalid e-w resolution field: 0

With OGR debugging on:

GRASS 7.3.svn (modis_sinusoidal):~ > PROJ_DEBUG=ON CPL_DEBUG=ON  v.import in=point_eu_laea.shp out=point_eu_laea2
GDAL: In GDALDestroy - unloading GDAL shared library.
Shape: DBF Codepage = LDID/87 for point_eu_laea.shp
Shape: Treating as encoding 'ISO-8859-1'.
GDAL: GDALOpen(point_eu_laea.shp, this=0x255f710) succeeds as ESRI Shapefile.
GDAL: In GDALDestroy - unloading GDAL shared library.
GDAL: force close of point_eu_laea.shp (0x255f710) in GDALDriverManager cleanup.
GDAL: GDALClose(point_eu_laea.shp, this=0x255f710)
Shape: DBF Codepage = LDID/87 for point_eu_laea.shp
Shape: Treating as encoding 'ISO-8859-1'.
GDAL: GDALOpen(point_eu_laea.shp, this=0x9d57e0) succeeds as ESRI Shapefile.
ERROR: Invalid e-w resolution field: 0
GDAL: In GDALDestroy - unloading GDAL shared library.
GDAL: force close of point_eu_laea.shp (0x9d57e0) in GDALDriverManager cleanup.
GDAL: GDALClose(point_eu_laea.shp, this=0x9d57e0)
ERROR: Unable to create location from OGR datasource <point_eu_laea.shp>

The projection of the location:

GRASS 7.3.svn (modis_sinusoidal):~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name       : unnamed
a          : 6371007.181
es         : 0
proj       : sinu
lon_0      : 0
x_0        : 0
y_0        : 0
no_defs    : defined
-PROJ_UNITS------------------------------------------------
unit       : Meter
units      : Meters
meters     : 1

I am a bit clueless here.

Attachments (4)

point_eu_laea.zip (1.0 KB ) - added by neteler 7 years ago.
Test point in EU LAEA (EPSG3035) projection - SHAPE file
v.in.ogr.diff (1018 bytes ) - added by neteler 7 years ago.
0.5 map units trick (dirty hack)
no_vector_default_region.diff (1.7 KB ) - added by mlennert 7 years ago.
diff for not basing region on vector extent
no_vector_default_region_2.diff (8.1 KB ) - added by mmetz 7 years ago.
fixing another diff for not basing region on vector extent

Download all attachments as: .zip

Change History (24)

by neteler, 7 years ago

Attachment: point_eu_laea.zip added

Test point in EU LAEA (EPSG3035) projection - SHAPE file

comment:1 by annakrat, 7 years ago

Possibly duplicate of #2878?

comment:2 by mlennert, 7 years ago

Interestingly:

g.proj georef=point_eu_laea.shp location=test2
Essai d'ouverture à l'aide d'OGR...
...réussi.
Location <test2> created
You can switch to the new location by
`g.mapset mapset=PERMANENT location=test2`

but

v.in.ogr point_eu_laea.shp location=test
ATTENTION: All available OGR layers will be imported into vector map
           <point_eu_laea>
Location <test> created
ERREUR :Invalid e-w resolution field: 0

Both are created, but:

diff -u test/PERMANENT/ test2/PERMANENT/
diff -u test/PERMANENT/.bashrc test2/PERMANENT/.bashrc
--- test/PERMANENT/.bashrc	2017-09-06 15:54:44.389963786 +0200
+++ test2/PERMANENT/.bashrc	2017-09-06 15:54:55.545841900 +0200
@@ -1,5 +1,5 @@
 test -r ~/.alias && . ~/.alias
-PS1='GRASS 7.3.svn (test):\w > '
+PS1='GRASS 7.3.svn (test2):\w > '
 grass_prompt() {
 	LOCATION="`g.gisenv get=GISDBASE,LOCATION_NAME,MAPSET separator='/'`"
 	if test -d "$LOCATION/grid3/G3D_MASK" && test -f "$LOCATION/cell/MASK" ; then
diff -u test/PERMANENT/DEFAULT_WIND test2/PERMANENT/DEFAULT_WIND
--- test/PERMANENT/DEFAULT_WIND	2017-09-06 15:51:39.891978100 +0200
+++ test2/PERMANENT/DEFAULT_WIND	2017-09-06 15:52:13.227614355 +0200
@@ -10,9 +10,9 @@
 n-s resol:  0
 top:        1.000000000000000
 bottom:     0.000000000000000
-cols3:      1
-rows3:      1
+cols3:      20
+rows3:      20
 depths:     1
-e-w resol3: 1
-n-s resol3: 1
+e-w resol3: 0
+n-s resol3: 0
 t-b resol:  1
Common subdirectories: test/PERMANENT/sqlite and test2/PERMANENT/sqlite
Only in test/PERMANENT/: .tmp
diff -u test/PERMANENT/WIND test2/PERMANENT/WIND
--- test/PERMANENT/WIND	2017-09-06 15:51:39.891978100 +0200
+++ test2/PERMANENT/WIND	2017-09-06 15:52:13.227614355 +0200
@@ -10,9 +10,9 @@
 n-s resol:  0
 top:        1.000000000000000
 bottom:     0.000000000000000
-cols3:      1
-rows3:      1
+cols3:      20
+rows3:      20
 depths:     1
-e-w resol3: 1
-n-s resol3: 1
+e-w resol3: 0
+n-s resol3: 0
 t-b resol:  1

and in both I get:

g.proj -p
ERREUR :Invalid e-w resolution field: 0

because DEFAULT_WIND/WIND contains:

proj:       99
zone:       0
north:      1205000
south:      1205000
east:       3674000
west:       3674000
cols:       20
rows:       20
e-w resol:  0
n-s resol:  0
top:        1.000000000000000
bottom:     0.000000000000000
cols3:      1
rows3:      1
depths:     1
e-w resol3: 1
n-s resol3: 1
t-b resol:  1

which seems logical as a point does not have dimensions, but the default region should still be defined with one cell and a resolution of 1....

So, yes, there seems to be a link to #2878...

in reply to:  2 ; comment:3 by mmetz, 7 years ago

Replying to mlennert:

[...] I get:

g.proj -p
ERREUR :Invalid e-w resolution field: 0

because DEFAULT_WIND/WIND contains:

proj:       99
zone:       0
north:      1205000
south:      1205000
east:       3674000
west:       3674000
[...]

which seems logical as a point does not have dimensions, but the default region should still be defined with one cell and a resolution of 1....

A simple solution would be to add/subtract 0.5 if north == south and/or east == west with rows and/or cols set to 1. For latlong, this is only guaranteed to work in trunk.

So, yes, there seems to be a link to #2878...

It seems so, because g.proj uses code copied from v.in.ogr.

in reply to:  3 ; comment:4 by mlennert, 7 years ago

Replying to mmetz:

Replying to mlennert:

[...] I get:

g.proj -p
ERREUR :Invalid e-w resolution field: 0

because DEFAULT_WIND/WIND contains:

proj:       99
zone:       0
north:      1205000
south:      1205000
east:       3674000
west:       3674000
[...]

which seems logical as a point does not have dimensions, but the default region should still be defined with one cell and a resolution of 1....

A simple solution would be to add/subtract 0.5 if north == south and/or east == west with rows and/or cols set to 1. For latlong, this is only guaranteed to work in trunk.

I think that the default region created when using a vector as input for georeference should always be the default n=1 s=0 e=1 w=0 res=1 region (i.e. the same that is created when creating a location without determining a default region). One could argue that the user would like to have a default region with the extents of the vector used for creating the location, but as there is no way to determine a sensible resolution, the risk is high that you get an arbitrarily large region with a resolution of 1, i.e. a very large number of pixels. Or we would have to decide on a reasonable default number of pixels and adjust the region accordingly. But I think that having an obviously non-sense region is better as it forces the user to explicitly set the region.

in reply to:  4 ; comment:5 by mmetz, 7 years ago

Replying to mlennert:

Replying to mmetz:

Replying to mlennert:

[...] I get:

g.proj -p
ERREUR :Invalid e-w resolution field: 0

because DEFAULT_WIND/WIND contains:

proj:       99
zone:       0
north:      1205000
south:      1205000
east:       3674000
west:       3674000
[...]

which seems logical as a point does not have dimensions, but the default region should still be defined with one cell and a resolution of 1....

A simple solution would be to add/subtract 0.5 if north == south and/or east == west with rows and/or cols set to 1. For latlong, this is only guaranteed to work in trunk.

I think that the default region created when using a vector as input for georeference should always be the default n=1 s=0 e=1 w=0 res=1 region (i.e. the same that is created when creating a location without determining a default region). One could argue that the user would like to have a default region with the extents of the vector used for creating the location, but as there is no way to determine a sensible resolution, the risk is high that you get an arbitrarily large region with a resolution of 1, i.e. a very large number of pixels.

I aggree.

Or we would have to decide on a reasonable default number of pixels and adjust the region accordingly.

I don't see a way to determine a reasonable resolution/number of pixels from vector extents.

But I think that having an obviously non-sense region is better as it forces the user to explicitly set the region.

Usually, users need to set (or at least review) the region anyway.

That means that duplicate code in g.proj/v.in.ogr where new region settings are estimated needs to be removed.

in reply to:  5 comment:6 by mlennert, 7 years ago

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to mlennert:

[...] I get:

g.proj -p
ERREUR :Invalid e-w resolution field: 0

because DEFAULT_WIND/WIND contains:

proj:       99
zone:       0
north:      1205000
south:      1205000
east:       3674000
west:       3674000
[...]

which seems logical as a point does not have dimensions, but the default region should still be defined with one cell and a resolution of 1....

A simple solution would be to add/subtract 0.5 if north == south and/or east == west with rows and/or cols set to 1. For latlong, this is only guaranteed to work in trunk.

I think that the default region created when using a vector as input for georeference should always be the default n=1 s=0 e=1 w=0 res=1 region (i.e. the same that is created when creating a location without determining a default region). One could argue that the user would like to have a default region with the extents of the vector used for creating the location, but as there is no way to determine a sensible resolution, the risk is high that you get an arbitrarily large region with a resolution of 1, i.e. a very large number of pixels.

I aggree.

Or we would have to decide on a reasonable default number of pixels and adjust the region accordingly.

I don't see a way to determine a reasonable resolution/number of pixels from vector extents.

But I think that having an obviously non-sense region is better as it forces the user to explicitly set the region.

Usually, users need to set (or at least review) the region anyway.

That means that duplicate code in g.proj/v.in.ogr where new region settings are estimated needs to be removed.

+1

by neteler, 7 years ago

Attachment: v.in.ogr.diff added

0.5 map units trick (dirty hack)

comment:7 by neteler, 7 years ago

I have attached the 0.5 suggestion being tested as dirty hack.

Seems to do the job while a correct solution in v.in.ogr and g.proj would obviously be preferred... (not sure how to implement it).

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

Replying to neteler:

I have attached the 0.5 suggestion being tested as dirty hack.

Seems to do the job while a correct solution in v.in.ogr and g.proj would obviously be preferred... (not sure how to implement it).

I don't know if you might have to make sure that cellhd.ns_res is big enough to avoid resolutions which might cause precision issues ? Also, in #2878 I mention issues with vertical resolution.

Personally, I think it better to go with the idea that when a location is defined based on a vector, no attempt is made to define the region based on the vector.

by mlennert, 7 years ago

diff for not basing region on vector extent

comment:9 by mlennert, 7 years ago

I just attached a diff which I think should work for not attempting to create a region based on vector extent. The region will simply be n=1 s=0 e=1 w=0 res=1.

I don't have time for much testing, so if others could, this would be good. Not sure though that this is secure enough to introduce now into the 7.2.2 release. Probably should be noted under known bugs and backported to release72 after the 7.2.2 release.

in reply to:  9 ; comment:10 by mmetz, 7 years ago

Replying to mlennert:

I just attached a diff which I think should work for not attempting to create a region based on vector extent. The region will simply be n=1 s=0 e=1 w=0 res=1.

For v.in.ogr, some more changes are needed, see attached patch no_vector_default_region_2.diff (the patch also solves some TODOs).

I don't have time for much testing, so if others could, this would be good. Not sure though that this is secure enough to introduce now into the 7.2.2 release. Probably should be noted under known bugs and backported to release72 after the 7.2.2 release.

+1

in reply to:  9 comment:11 by neteler, 7 years ago

Replying to mlennert:

I don't have time for much testing, so if others could, this would be good. Not sure though that this is secure enough to introduce now into the 7.2.2 release. Probably should be noted under known bugs and backported to release72 after the 7.2.2 release.

... this means that we will have to release 7.2.3 for sure because that's a relevant bug.

in reply to:  10 ; comment:12 by neteler, 7 years ago

Replying to mmetz:

For v.in.ogr, some more changes are needed, see attached patch no_vector_default_region_2.diff (the patch also solves some TODOs).

I get a compilation error:

[mneteler@oboe v.in.ogr ]$ make
gcc  -O2 -march=native -std=gnu99 -fexceptions -fstack-protector -m64 -fdiagnostics-color   -I/home/mneteler/software/grass73/dist.x86_64-pc-linux-gnu/include -I/home/mneteler/software/grass73/dist.x86_64-pc-linux-gnu/include   -I/usr/include/gdal -I/usr/include -DPACKAGE=\""grassmods"\"   -I/usr/include/gdal -I/home/mneteler/software/grass73/dist.x86_64-pc-linux-gnu/include -I/home/mneteler/software/grass73/dist.x86_64-pc-linux-gnu/include -DRELDIR=\"vector/v.in.ogr\" -o OBJ.x86_64-pc-linux-gnu/main.o -c main.c
main.c: In function ‘main’:
main.c:886:35: error: ‘have_ogr_extents’ undeclared (first use in this function); did you mean ‘have_ogr_extent’?
     if (flag.no_clean->answer || !have_ogr_extents) {
                                   ^~~~~~~~~~~~~~~~
                                   have_ogr_extent
main.c:886:35: note: each undeclared identifier is reported only once for each function it appears in
make: *** [../../include/Make/Compile.make:32: OBJ.x86_64-pc-linux-gnu/main.o] Error 1

by mmetz, 7 years ago

fixing another diff for not basing region on vector extent

in reply to:  12 comment:13 by mmetz, 7 years ago

Replying to neteler:

Replying to mmetz:

For v.in.ogr, some more changes are needed, see attached patch no_vector_default_region_2.diff (the patch also solves some TODOs).

I get a compilation error:

Oops, fixed in the updated attached patch no_vector_default_region_2.diff.

About backporting: v.in.ogr in trunk is already quite different from relbr72, and with this patch, it would be more different. IMHO, the fix for OSM import and this fix are both bug fixes, i.e. once confirmed working, patched v.in.ogr in trunk should be backported as is to relbr72.

comment:14 by neteler, 7 years ago

Thanks for the new updated no_vector_default_region_2.diff, tested here:

GRASS 7.3.svn (modis_sinusoidal):~ > v.import --v in=~/point_eu_laea.shp out=point_eu_laea2Creating temporary location for </home/mundialis/point_eu_laea.shp>...
Using OGR driver 'ESRI Shapefile'
Projection of dataset does not appear to match current location.

GRASS LOCATION PROJ_INFO is:
name: unnamed
a: 6371007.181
es: 0
proj: sinu
lon_0: 0
x_0: 0
y_0: 0
no_defs: defined

Import dataset PROJ_INFO is:
name: ETRS89_ETRS_LAEA
datum: etrs89
ellps: grs80
proj: laea
lat_0: 52
lon_0: 10
x_0: 4321000
y_0: 3210000
no_defs: defined

In case of no significant differences in the projection definitions, use
the -o flag to ignore them and use current location definition.
Consider generating a new location with 'location' parameter from input
data set.
switch to temp location
-PROJ_INFO-------------------------------------------------
name : ETRS89_ETRS_LAEA
datum : etrs89
ellps : grs80
proj : laea
lat_0 : 52
lon_0 : 10
x_0 : 4321000
y_0 : 3210000
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : Meter
units : Meters
meters : 1
Importing </home/mundialis/point_eu_laea.shp> ...
Using OGR driver 'ESRI Shapefile'
Projection of input dataset and current location appear to match
Check if OGR layer <point_eu_laea> contains polygons...
 100%
Importing 1 features
Using native format
Default driver / database set to:
driver: sqlite
database: $GISDBASE/$LOCATION_NAME/$MAPSET/sqlite/sqlite.db
Column name <cat> renamed to <cat_>
Importing 1 features (OGR layer <point_eu_laea>)...
 100%
-----------------------------------------------------
Building topology for vector map <point_eu_laea2@PERMANENT>...
Registering primitives...
One primitive registered
One vertex registered
Building areas...
 100%
0 areas built
0 isles built
Attaching islands...
Attaching centroids...
 100%
Topology was built
Number of nodes: 0
Number of primitives: 1
Number of points: 1
Number of lines: 0
Number of boundaries: 0
Number of centroids: 0
Number of areas: 0
Number of isles: 0
Reprojecting <point_eu_laea2>...
Input Projection Parameters:  +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +no_defs +a=6378137 +rf=298.257222101 +towgs84=0.000,0.000,0.000
Input Unit Factor: 1
Output Projection Parameters:  +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +no_defs +a=6371007.181 +b=6371007.181
Output Unit Factor: 1
Using native format
Reprojecting primitives ...
Building topology for vector map
<point_eu_laea2@Germany_LST_1km_daily_gdd_2016>...
Registering primitives...
One primitive registered
One vertex registered
Building areas...
 100%
0 areas built
0 isles built
Attaching islands...
Attaching centroids...
 100%
Topology was built
Number of nodes: 0
Number of primitives: 1
Number of points: 1
Number of lines: 0
Number of boundaries: 0
Number of centroids: 0
Number of areas: 0
Number of isles: 0

Now the file is smoothly imported :-)

comment:15 by neteler, 7 years ago

Markus M, do you plan to submit the patch?

in reply to:  15 ; comment:16 by mmetz, 7 years ago

Replying to neteler:

Markus M, do you plan to submit the patch?

Yes, I plan to submit the patch, but the patch affects the options spatial=, where= and the -r and -e flags which I want to test a bit more before submitting.

in reply to:  16 comment:17 by mmetz, 7 years ago

Replying to mmetz:

Replying to neteler:

Markus M, do you plan to submit the patch?

Yes, I plan to submit the patch, but the patch affects the options spatial=, where= and the -r and -e flags which I want to test a bit more before submitting.

Fixed in trunk with r71490,1 (v.in.ogr, g.proj).

comment:18 by neteler, 7 years ago

Milestone: 7.2.27.2.3

Ticket retargeted after milestone closed

comment:19 by martinl, 6 years ago

What is status of this issue? Backport needed?

in reply to:  19 comment:20 by neteler, 6 years ago

Milestone: 7.2.37.4.0
Resolution: fixed
Status: newclosed
Version: svn-trunksvn-releasebranch74

Replying to martinl:

What is status of this issue? Backport needed?

Since it was fixed in trunk 6 months ago it should be part of 7.4.0, that's fine. Closing.

Note: See TracTickets for help on using tickets.