r.in.gdal/v.in.ogr: reprojection at import

[mlennert]

I would like to float the idea, in order to kick off a discussion: how difficult would it be to implement reprojection at import, both in r.in.gdal and v.in.ogr.

The simplest case would be the import of a data source with known and gdal/ogr-recognizable projection system into a location with a different projection. If the data source does not contain the info about its projection, then the user would have to add that info.

This should complement, not replace the current system of reprojecting from one location to another as this latter system allows more subtle fine-tuning of the process.

Moritz

[martinl]

I was thinking about similar issue relating to wxGUI data catalog. The interactive tool which could enable user to rename, remove or copy maps from different mapsets and also locations (this step would contain also reprojection in the case that locations have different projections) or import data in external format (GDAL/OGR based) including a preview. There is a prototype of such application (1) but it will be probably never integrated into wxGUI core, in other words must be implemented from the scratch.

Back to the topic, reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

(1) ​http://grasswiki.osgeo.org/wiki/WxPython-based_GUI_for_GRASS#GRASS_Catalog

[mmetz]

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

Markus M

[mlennert]

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

Moritz

[sbl]

See also: #2486

[mmetz]

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

For example, it is mathematically not possible to reproject a global map in latlong to UTM or most other true projection systems. You should only reproject a subregion (that part that corresponds to the current region in the target projection) to the target projection. That also for raster data that the current region first needs to be adjusted according to the desired resolution. Not to mention the resampling method to be used for reprojection. The default is nearest neighbor, but nn usually implies information loss. Surfaces are best reprojected with bilinear or bicubic, imagery best with e.g. lanczos. This can not be determined automatically. For high-precision reprojection, the target resolution would depend not only on the source resolution, but also on the resampling method (think e.g. MODIS quality assessment + pixel shift after reprojection).

For vector data, lines and boundaries need to be densified first (adding additional vertices). I have implemented this for v.proj in G7. Still, the reprojection should be limited to the target current region (currently possible with v.in.region + v.proj + v.overlay + v.proj) because it is mathematically not possible to reproject e.g. a global map in latlong to UTM or most other true projection systems.

[mmetz]

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource. It is by purpose educational, i.e. the user needs to specify the resampling method and by default the average target resolution is estimated, but the input is not reprojected. This forces the user to think about both the resampling method and a reasonable output resolution. If in doubt, read the manual...

[annakrat]

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource. It is by purpose educational, i.e. the user needs to specify the resampling method and by default the average target resolution is estimated, but the input is not reprojected. This forces the user to think about both the resampling method and a reasonable output resolution. If in doubt, read the manual...

Thanks, that's great! Eventually, this should be a core module and should be incorporated in the import dialog. Although the gui should be simple, we can still force users to actively decide the method and resolution.

[mlennert]

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource. It is by purpose educational, i.e. the user needs to specify the resampling method and by default the average target resolution is estimated, but the input is not reprojected. This forces the user to think about both the resampling method and a reasonable output resolution. If in doubt, read the manual...

Great work, thanks a lot (including for the educational approach) !

A few little bugs I've come upon:

1) 'output' is said to be =input by default but that does not seem to be the case anywhere in the source code, so output is a mandatory option.

2) There are some strings (n,s,e,w) that need to be converted to float for calculation in line ​290.

3) The way I would understand the 'extents' option this should import the entire map. So this would entail that the region is adapted to the map prior to reprojection. The '-n' flag of r.proj does not do this, IIUC, so when I try to import I get an empty map. I guess this needs a first dry run of r.proj with the -g flag.

Thanks again for this !

Moritz

[mmetz]

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource. It is by purpose educational, i.e. the user needs to specify the resampling method and by default the average target resolution is estimated, but the input is not reprojected. This forces the user to think about both the resampling method and a reasonable output resolution. If in doubt, read the manual...

Great work, thanks a lot (including for the educational approach) !

A few little bugs I've come upon:

1) 'output' is said to be =input by default but that does not seem to be the case anywhere in the source code, so output is a mandatory option.

2) There are some strings (n,s,e,w) that need to be converted to float for calculation in line ​290.

3) The way I would understand the 'extents' option this should import the entire map. So this would entail that the region is adapted to the map prior to reprojection. The '-n' flag of r.proj does not do this, IIUC, so when I try to import I get an empty map. I guess this needs a first dry run of r.proj with the -g flag.

Thanks for testing. All fixed in r64352.

[mmetz]

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource.

Now there is also v.in.proj which imports and reprojects (if necessary) a OGR datasource. It is still a bit noisy (lots of messages if direct import fails).

[mlennert]

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource. It is by purpose educational, i.e. the user needs to specify the resampling method and by default the average target resolution is estimated, but the input is not reprojected. This forces the user to think about both the resampling method and a reasonable output resolution. If in doubt, read the manual...

Great work, thanks a lot (including for the educational approach) !

A few little bugs I've come upon:

1) 'output' is said to be =input by default but that does not seem to be the case anywhere in the source code, so output is a mandatory option.

2) There are some strings (n,s,e,w) that need to be converted to float for calculation in line ​290.

3) The way I would understand the 'extents' option this should import the entire map. So this would entail that the region is adapted to the map prior to reprojection. The '-n' flag of r.proj does not do this, IIUC, so when I try to import I get an empty map. I guess this needs a first dry run of r.proj with the -g flag.

Thanks for testing. All fixed in r64352.

I still get an error:

Traceback (most recent call last):
  File
"/data/home/mlennert/SRC/GRASS/grass_trunk/dist.x86_64
-unknown-linux-gnu/scripts/r.in.proj", line 337, in <module>
    main()
  File
"/data/home/mlennert/SRC/GRASS/grass_trunk/dist.x86_64
-unknown-linux-gnu/scripts/r.in.proj", line 302, in main
    grass.message(_("Specified target resolution: %g") %
outres)
TypeError: float argument required, not str

I can solve this with

===================================================================
--- r.in.proj.py	(révision 64355)
+++ r.in.proj.py	(copie de travail)
@@ -231,7 +231,7 @@
     region = grass.region()
     
     if tgtres is not None:
-        outres = tgtres
+        outres = float(tgtres)
     else:
         outres = (region['ewres'] + region['nsres']) / 2.0

but don't know if that is your intention.

Two more usability remarks:

1. You actually run r.in.gdal (line 187) and you test afterwards, whether a projection system was recognized or an XY-location created. Could this test be put after line 160 by looking at 'insrs' ? This would avoid going through the import which can take a while if the dataset is heavy.

2. When the '-i' flag is not set a short message at the end should inform the user about the fact that the layer has not been imported, especially since the module spits out a 'Raster map <XYZ> created' during the process (coming from the creation of the raster map in the temp loc) and that you have to dig into the doc to understand that you have to explicitly set the -i flag.

[mlennert]

Replying to mmetz:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource.

Now there is also v.in.proj which imports and reprojects (if necessary) a OGR datasource. It is still a bit noisy (lots of messages if direct import fails).

Great, but you put it into the raster directory. Shouldn't it be in the vector directory ?

[mlennert]

Replying to mmetz:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource.

Now there is also v.in.proj which imports and reprojects (if necessary) a OGR datasource. It is still a bit noisy (lots of messages if direct import fails).

Great. Seems to work without problems in simple test.

One question: why did you decide to not use the g.proj -j comparison test before direct import here ? The current error message from the failed direct import is a bit disconcerting.

[mmetz]

Replying to mlennert:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource. It is by purpose educational, i.e. the user needs to specify the resampling method and by default the average target resolution is estimated, but the input is not reprojected. This forces the user to think about both the resampling method and a reasonable output resolution. If in doubt, read the manual...

Great work, thanks a lot (including for the educational approach) !

A few little bugs I've come upon:

1) 'output' is said to be =input by default but that does not seem to be the case anywhere in the source code, so output is a mandatory option.

2) There are some strings (n,s,e,w) that need to be converted to float for calculation in line ​290.

3) The way I would understand the 'extents' option this should import the entire map. So this would entail that the region is adapted to the map prior to reprojection. The '-n' flag of r.proj does not do this, IIUC, so when I try to import I get an empty map. I guess this needs a first dry run of r.proj with the -g flag.

Thanks for testing. All fixed in r64352.

I still get an error:

[...]

Fixed.

Two more usability remarks:

1. You actually run r.in.gdal (line 187) and you test afterwards, whether a projection system was recognized or an XY-location created. Could this test be put after line 160 by looking at 'insrs' ? This would avoid going through the import which can take a while if the dataset is heavy.

OK, done in r64358. First come the CRS tests, then the actual import.

2. When the '-i' flag is not set a short message at the end should inform the user about the fact that the layer has not been imported, especially since the module spits out a 'Raster map <XYZ> created' during the process (coming from the creation of the raster map in the temp loc) and that you have to dig into the doc to understand that you have to explicitly set the -i flag.

The module now says "The input <XYZ> can be imported and reprojected with the -i flag".

[mmetz]

Replying to mlennert:

Replying to mmetz:

Replying to mmetz:

Replying to mmetz:

Replying to mlennert:

Replying to mmetz:

Replying to martinl:

reprojection option available for r.in.gdal/v.in.ogr would a step forward from user perspective.

You would need to supply all options of [r|v].proj. Moreover, vector data might need preprocessing (adding vertices to lines/boundaries) in order to provide realistic and topologically correct results. The reprojection option could IMHO be best implemented in a script which calls r.in.gdal/v.in.ogr and then r.proj/(v.split + v.proj).

I had actually thought that it might be possible to integrate GDAL/OGR library tools such as GDALWarpOperation and OGRCoordinateTransformation directly into r.in.gdal/v.in.ogr.

But maybe you're right and we should not touch the general GRASS structure of import into one location + reproject into another, and rather use a wrapper.

(Ongoing discussion on the dev-ml)

My point for a wrapper script which calls r.in.gdal/v.in.ogr and then r.proj/v.proj is about reprojection pitfalls.

I have created the addon script r.in.proj which imports and reprojects (if necessary) a GDAL datasource.

Now there is also v.in.proj which imports and reprojects (if necessary) a OGR datasource. It is still a bit noisy (lots of messages if direct import fails).

Great. Seems to work without problems in simple test.

One question: why did you decide to not use the g.proj -j comparison test before direct import here ? The current error message from the failed direct import is a bit disconcerting.

I have changed this and synced r.in.proj with v.in.proj in r64359. First, a temporary location is created from the input, this is fast. Then the projections are compared (g.proj -j output) and, if equal, the input is directly imported without reprojection. Otherwise, the input is imported in the just created temporary location and from there reprojected to the current location.

[martinl]

Replying to annakrat:

Thanks, that's great! Eventually, this should be a core module and should be incorporated in the import dialog. Although the gui should be simple, we can still force users to actively decide the method and resolution.

I absolutely agree, are you planning to do that? Thanks, Martin

[mmetz]

Replying to martinl:

Replying to annakrat:

Thanks, that's great! Eventually, this should be a core module and should be incorporated in the import dialog. Although the gui should be simple, we can still force users to actively decide the method and resolution.

I absolutely agree, are you planning to do that? Thanks, Martin

There is a GUI-related issue. The two modules currently have each two input options, input_file and input_directory. It would be nice if these two options could be combined in one, key input, type string, and the GUI would provide an interface similar to r.in.gdal/v.in.ogr where you can select if you want to import a file or a directory or something from a database. The difference to the current r.in.gdal/v.in.ogr interface would be that input is only one GDAL/OGR datasource. Each GDAL/OGR datasource may contain several band/layers, but that is handled by r.in.gdal/v.in.ogr. GDAL supports currently 135 different formats, OGR 80 different formats. Not all of them are file based, thus the need for flexible specification of the input datasource.

The aim is to provide a simple yet powerful interface for users to import data which are (semi-)automatically reprojected if need be.

[martinl]

There are new modules r.import and v.import. Closing.