[GRASS-user] Re-projected Data Position Mismatch

GRASS GRASS user
,
1

Greetings,

I am attempting to merge data in a couple different projections into a single projection. But as you can see in the following image, the vector data does not align with the raster tiles. The horizontal positions match, but they don’t match vertically.

http://3DTOPO.com/misalignedData.jpg

I’ve never ran into this issue before using GRASS after merging countless projections together.

All of the data shown in the image is in “epsg: 3857” projection, here is the contents of its PROJ_INFO file:

name: WGS 84 / Pseudo-Mercator
datum: wgs84
ellps: wgs84
proj: merc
lat_ts: 0.0
lon_0: 0.0
x_0: 0.0
y_0: 0
k: 1.0
wktext: defined
no_defs: defined
towgs84: 0.000,0.000,0.000

I will refer to this location as the Comp Location.

If I import a raster tile as a new location, it has this as a resulting PROJ_INFO file:

name: WGS 84 / Pseudo-Mercator
datum: wgs84
ellps: wgs84
proj: merc
lat_ts: 0.0
lon_0: 0.0
x_0: 0.0
y_0: 0
k: 1.0
nadgrids: @null
wktext: defined
no_defs: defined

The only difference that I can tell is that the raster tile location has a "nadgrids: @null” line which I had to delete to import the raster tiles into the first location with r.proj, the Comp Location has a towgs84 transformation (but is 0,0,0), and no PROJ_EPSG file exists in the Comp Location. Alternatively I can import the raster tiles directly into the Comp Location using r.in.gdal since the projections match (I don’t get any mismatched projection warnings). In both cases, the raster tiles are located in the same position.

The PROJ_INFO file from the locations that the vector orange grid and green bounds was imported from is as follows:

name: GCS_North_American_1983
datum: nad83
ellps: grs80
proj: ll
no_defs: defined

I had imported the orange grid and green bounds from two different locations so I could try to determine what data was in the wrong position (raster tiles or vector) but, but the both the vector layers imported are in the same lonlat projection, so I guess is wasn’t a really good test, thus I am not entirely sure which data is in the correct location.

Any suggestions would be greatly appreciated.

Thanks,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

2

On Aug 22, 2017, at 12:16 PM, Jeshua Lacock <jeshua@3DTOPO.com> wrote:

and no PROJ_EPSG file exists in the Comp Location

Correction; no PROJ_EPSG file exists in the rater tile location, it does exist in the Comp Location.

3

* Jeshua Lacock <jeshua@3DTOPO.com> [2017-08-22 12:16:56 -0600]:

Greetings,

I am attempting to merge data in a couple different projections into a
single projection. But as you can see in the following image, the
vector data does not align with the raster tiles. The horizontal
positions match, but they don’t match vertically.

http://3DTOPO.com/misalignedData.jpg

Hello Jeshua,

the vector map seems to be a regular grid. If you want it to fit
exactly, can't you recreate it? Using `v.mkgrid` and cut it based on the
raster map(s) of interest?

Nikos

I’ve never ran into this issue before using GRASS after merging countless projections together.

All of the data shown in the image is in “epsg: 3857” projection, here is the contents of its PROJ_INFO file:

name: WGS 84 / Pseudo-Mercator
datum: wgs84
ellps: wgs84
proj: merc
lat_ts: 0.0
lon_0: 0.0
x_0: 0.0
y_0: 0
k: 1.0
wktext: defined
no_defs: defined
towgs84: 0.000,0.000,0.000

I will refer to this location as the Comp Location.

If I import a raster tile as a new location, it has this as a resulting PROJ_INFO file:

name: WGS 84 / Pseudo-Mercator
datum: wgs84
ellps: wgs84
proj: merc
lat_ts: 0.0
lon_0: 0.0
x_0: 0.0
y_0: 0
k: 1.0
nadgrids: @null
wktext: defined
no_defs: defined

The only difference that I can tell is that the raster tile location has a "nadgrids: @null” line which I had to delete to import the raster tiles into the first location with r.proj, the Comp Location has a towgs84 transformation (but is 0,0,0), and no PROJ_EPSG file exists in the Comp Location. Alternatively I can import the raster tiles directly into the Comp Location using r.in.gdal since the projections match (I don’t get any mismatched projection warnings). In both cases, the raster tiles are located in the same position.

The PROJ_INFO file from the locations that the vector orange grid and green bounds was imported from is as follows:

name: GCS_North_American_1983
datum: nad83
ellps: grs80
proj: ll
no_defs: defined

I had imported the orange grid and green bounds from two different locations so I could try to determine what data was in the wrong position (raster tiles or vector) but, but the both the vector layers imported are in the same lonlat projection, so I guess is wasn’t a really good test, thus I am not entirely sure which data is in the correct location.

Any suggestions would be greatly appreciated.

Thanks,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

--
Nikos Alexandris | Remote Sensing & Geomatics
GPG Key Fingerprint 6F9D4506F3CA28380974D31A9053534B693C4FB3

4

On Aug 22, 2017, at 1:03 PM, Nikos Alexandris <nik@nikosalexandris.net> wrote:

the vector map seems to be a regular grid. If you want it to fit
exactly, can't you recreate it? Using `v.mkgrid` and cut it based on the
raster map(s) of interest?

Hi Nikos,

Thanks for your suggestion. It is a regular grid, but I need to look up the names of some of the rectangles (which are attributes of the grid). Also, at some point I have to get all the vector data I have to match up with the raster tiles…

Best,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

5

To hopefully help troubleshoot; I just reprojected one of the raster tiles (from epsg: 3857), into the source location of one of the lonlat vectors (reverse projections from my OP), and the datasets are offset by the same distances. Since the dimensions of the raster are being changed (by r.proj), it leads me to think it must be a datum or coordinate system misalignment and not a projection issue.

For example, here is a test raster tile in the Comp Projection (epsg: 3857): http://3DTOPO.com/epsg_3857.png

And in the reverse lonlat projection: http://3DTOPO.com/lonlat.png

(the raster tile should be located in the red area in both examples)

6

On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com> wrote:

To hopefully help troubleshoot; I just reprojected one of the raster tiles (from epsg: 3857), into the source location of one of the lonlat vectors (reverse projections from my OP), and the datasets are offset by the same distances. Since the dimensions of the raster are being changed (by r.proj), it leads me to think it must be a datum or coordinate system misalignment and not a projection issue.

I have the same problem, working with NAIP imagery. It is related to:
https://trac.osgeo.org/grass/ticket/2229

I have to remove nadgrids: @null from the PROJ_INFO file to be able to
reproject into that location, but then it is shifted. gdalwarp helps.

For example, here is a test raster tile in the Comp Projection (epsg: 3857): http://3DTOPO.com/epsg_3857.png

And in the reverse lonlat projection: http://3DTOPO.com/lonlat.png

(the raster tile should be located in the red area in both examples)
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

7

On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna@gmail.com> wrote:

On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com> wrote:

To hopefully help troubleshoot; I just reprojected one of the raster tiles (from epsg: 3857), into the source location of one of the lonlat vectors (reverse projections from my OP), and the datasets are offset by the same distances. Since the dimensions of the raster are being changed (by r.proj), it leads me to think it must be a datum or coordinate system misalignment and not a projection issue.

I have the same problem, working with NAIP imagery. It is related to:
https://trac.osgeo.org/grass/ticket/2229

I have to remove nadgrids: @null from the PROJ_INFO file to be able to
reproject into that location, but then it is shifted. gdalwarp helps.

Hi Anna!

Thank you very much for confirming that! I am working with the NAIP imagery as well. :slight_smile:

I have found that their original Geotiff assets work perfectly.

In fact, I was very happy to stumble on to the fact that the complete NAIP archive (~250 terabytes) is available as a bucket drive on Amazon Web Services (AWS). So I setup GRASS instances to process the tiles, then download the processed, reprojected images compressed as JP2s. I am paying for the bandwidth and compute time, but I think its worth it for my purposes. I’ll be able to process and download the imagery I need in about 60 days compared to over 300 days without AWS!

Cheers,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

8

On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua@3dtopo.com> wrote:

On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna@gmail.com> wrote:

On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com> wrote:

To hopefully help troubleshoot; I just reprojected one of the raster tiles (from epsg: 3857), into the source location of one of the lonlat vectors (reverse projections from my OP), and the datasets are offset by the same distances. Since the dimensions of the raster are being changed (by r.proj), it leads me to think it must be a datum or coordinate system misalignment and not a projection issue.

I have the same problem, working with NAIP imagery. It is related to:
https://trac.osgeo.org/grass/ticket/2229

I have to remove nadgrids: @null from the PROJ_INFO file to be able to
reproject into that location, but then it is shifted. gdalwarp helps.

EPSG:3857 is problematic because it
“Uses spherical development of ellipsoidal coordinates. Relative to WGS 84 / World Mercator (CRS code 3395) errors of 0.7 percent in scale and differences in northing of up to 43km in the map (equivalent to 21km on the ground) may arise.”

Therefore you would need to use the WKT EXTENSION PROJ4:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +wktext +no_defs

without +nadgrids=@null

In GRASS, you could use this proj4 definition to create a new location, then import the data (the -o flag might be needed), then reproject to the desired CRS.

Considering this particular CRS (EPSG:3857), it is strange than gdalwarp works while GRASS reprojection does not work because GRASS uses GDAL to convert WKT to proj4, then to GRASS terminology. Some debugging is needed to figure out why within GRASS, the conversion from WKT to proj4 (using GDAL) produces wrong results.

Markus M

Hi Anna!

Thank you very much for confirming that! I am working with the NAIP imagery as well. :slight_smile:

I have found that their original Geotiff assets work perfectly.

In fact, I was very happy to stumble on to the fact that the complete NAIP archive (~250 terabytes) is available as a bucket drive on Amazon Web Services (AWS). So I setup GRASS instances to process the tiles, then download the processed, reprojected images compressed as JP2s. I am paying for the bandwidth and compute time, but I think its worth it for my purposes. I’ll be able to process and download the imagery I need in about 60 days compared to over 300 days without AWS!

Cheers,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

9

On Wed, Aug 30, 2017 at 10:46 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:
...

Considering this particular CRS (EPSG:3857), it is strange than gdalwarp
works while GRASS reprojection does not work because GRASS uses GDAL to
convert WKT to proj4, then to GRASS terminology. Some debugging is needed to
figure out why within GRASS, the conversion from WKT to proj4 (using GDAL)
produces wrong results.

For whom is interested to do so, here explanations on how to debug all
this (we crafted this part of the README at the FOSS4G 2016 in Bonn
during the code sprint):
https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/README.txt

markusN

10

On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua@3dtopo.com> wrote:

> On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna@gmail.com>
> wrote:
>
> On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com>
> wrote:
>>
>> To hopefully help troubleshoot; I just reprojected one of the raster
>> tiles (from epsg: 3857), into the source location of one of the lonlat
>> vectors (reverse projections from my OP), and the datasets are offset by the
>> same distances. Since the dimensions of the raster are being changed (by
>> r.proj), it leads me to think it must be a datum or coordinate system
>> misalignment and not a projection issue.
>
> I have the same problem, working with NAIP imagery. It is related to:
> https://trac.osgeo.org/grass/ticket/2229
>
> I have to remove nadgrids: @null from the PROJ_INFO file to be able to
> reproject into that location, but then it is shifted. gdalwarp helps.

EPSG:3857 is problematic because it
"Uses spherical development of ellipsoidal coordinates. Relative to WGS 84 /
World Mercator (CRS code 3395) errors of 0.7 percent in scale and
differences in northing of up to 43km in the map (equivalent to 21km on the
ground) may arise."

Therefore you would need to use the WKT EXTENSION PROJ4:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +wktext +no_defs

without +nadgrids=@null

In GRASS, you could use this proj4 definition to create a new location, then
import the data (the -o flag might be needed), then reproject to the desired
CRS.

Considering this particular CRS (EPSG:3857), it is strange than gdalwarp
works while GRASS reprojection does not work because GRASS uses GDAL to
convert WKT to proj4, then to GRASS terminology. Some debugging is needed to
figure out why within GRASS, the conversion from WKT to proj4 (using GDAL)
produces wrong results.

Thank you, yes that works. I looked at lib/proj/convert.c where I
assume the problem is. If I interpret it correctly it basically
discards +a and +b from the EXTENSION and instead picks WGS84
ellipsoid because of wkt

DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]]

But I don't really know what that extension actually means or how
GRASS should treat it.

This is the WKT which is processed:

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0],EXTENSION["PROJ4","+proj=merc
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +wktext +no_defs"]]

Anna

Markus M

Hi Anna!

Thank you very much for confirming that! I am working with the NAIP
imagery as well. :slight_smile:

I have found that their original Geotiff assets work perfectly.

In fact, I was very happy to stumble on to the fact that the complete NAIP
archive (~250 terabytes) is available as a bucket drive on Amazon Web
Services (AWS). So I setup GRASS instances to process the tiles, then
download the processed, reprojected images compressed as JP2s. I am paying
for the bandwidth and compute time, but I think its worth it for my
purposes. I’ll be able to process and download the imagery I need in about
60 days compared to over 300 days without AWS!

Cheers,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

11

On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová <kratochanna@gmail.com> wrote:

On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua@3dtopo.com> wrote:

On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna@gmail.com>
wrote:

On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com>
wrote:

To hopefully help troubleshoot; I just reprojected one of the raster
tiles (from epsg: 3857), into the source location of one of the lonlat
vectors (reverse projections from my OP), and the datasets are offset by the
same distances. Since the dimensions of the raster are being changed (by
r.proj), it leads me to think it must be a datum or coordinate system
misalignment and not a projection issue.

I have the same problem, working with NAIP imagery. It is related to:
https://trac.osgeo.org/grass/ticket/2229

I have to remove nadgrids: @null from the PROJ_INFO file to be able to
reproject into that location, but then it is shifted. gdalwarp helps.

EPSG:3857 is problematic because it
“Uses spherical development of ellipsoidal coordinates. Relative to WGS 84 /
World Mercator (CRS code 3395) errors of 0.7 percent in scale and
differences in northing of up to 43km in the map (equivalent to 21km on the
ground) may arise.”

Therefore you would need to use the WKT EXTENSION PROJ4:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +wktext +no_defs

without +nadgrids=@null

In GRASS, you could use this proj4 definition to create a new location, then
import the data (the -o flag might be needed), then reproject to the desired
CRS.

Considering this particular CRS (EPSG:3857), it is strange than gdalwarp
works while GRASS reprojection does not work because GRASS uses GDAL to
convert WKT to proj4, then to GRASS terminology. Some debugging is needed to
figure out why within GRASS, the conversion from WKT to proj4 (using GDAL)
produces wrong results.

Thank you, yes that works. I looked at lib/proj/convert.c where I
assume the problem is. If I interpret it correctly it basically
discards +a and +b from the EXTENSION and instead picks WGS84
ellipsoid because of wkt

DATUM[“D_WGS_1984”,SPHEROID[“WGS_1984”,6378137.0,298.257223563]]

But I don’t really know what that extension actually means or how
GRASS should treat it.

This is the WKT which is processed:

PROJCS[“WGS_1984_Web_Mercator_Auxiliary_Sphere”,GEOGCS[“GCS_WGS_1984”,DATUM[“D_WGS_1984”,SPHEROID[“WGS_1984”,6378137.0,298.257223563]],PRIMEM[“Greenwich”,0.0],UNIT[“Degree”,0.017453292519943295]],PROJECTION[“Mercator_Auxiliary_Sphere”],PARAMETER[“False_Easting”,0.0],PARAMETER[“False_Northing”,0.0],PARAMETER[“Central_Meridian”,0.0],PARAMETER[“Standard_Parallel_1”,0.0],PARAMETER[“Auxiliary_Sphere_Type”,0.0],UNIT[“Meter”,1.0],EXTENSION[“PROJ4”,“+proj=merc
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +wktext +no_defs”]]

This WKT is processed correctly, GDAL converts this with OSRExportToProj4() to

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs

the only problem is +nadgrids=@null.

Unfortunately, GPJ_osr_to_grass() discards any datum and/or ellipsoid info in the proj4 term at
https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/convert.c#L477

Datum and/or ellipsoid definitions are determined later on from the original OGR spatial reference, which causes problems for EPSG:3857 and possibly other CRS’s.

The reason for this special handling is that GRASS has its own datum and ellipsoid definitions in lib/gis, i.e. in datum.table, datumtransform.table, ellipse.table, ellipse.table.solar.system.

It seems like a review of GRASS handling of GDAL CRS definitions is needed…

Markus M

Anna

Markus M

Hi Anna!

Thank you very much for confirming that! I am working with the NAIP
imagery as well. :slight_smile:

I have found that their original Geotiff assets work perfectly.

In fact, I was very happy to stumble on to the fact that the complete NAIP
archive (~250 terabytes) is available as a bucket drive on Amazon Web
Services (AWS). So I setup GRASS instances to process the tiles, then
download the processed, reprojected images compressed as JP2s. I am paying
for the bandwidth and compute time, but I think its worth it for my
purposes. I’ll be able to process and download the imagery I need in about
60 days compared to over 300 days without AWS!

Cheers,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

12

On Thu, Aug 31, 2017 at 3:47 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová <kratochanna@gmail.com>
wrote:

On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:
>
>
> On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua@3dtopo.com>
> wrote:
>>
>>
>> > On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna@gmail.com>
>> > wrote:
>> >
>> > On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com>
>> > wrote:
>> >>
>> >> To hopefully help troubleshoot; I just reprojected one of the raster
>> >> tiles (from epsg: 3857), into the source location of one of the
>> >> lonlat
>> >> vectors (reverse projections from my OP), and the datasets are
>> >> offset by the
>> >> same distances. Since the dimensions of the raster are being changed
>> >> (by
>> >> r.proj), it leads me to think it must be a datum or coordinate
>> >> system
>> >> misalignment and not a projection issue.
>> >
>> > I have the same problem, working with NAIP imagery. It is related to:
>> > https://trac.osgeo.org/grass/ticket/2229
>> >
>> > I have to remove nadgrids: @null from the PROJ_INFO file to be able
>> > to
>> > reproject into that location, but then it is shifted. gdalwarp helps.
>
> EPSG:3857 is problematic because it
> "Uses spherical development of ellipsoidal coordinates. Relative to WGS
> 84 /
> World Mercator (CRS code 3395) errors of 0.7 percent in scale and
> differences in northing of up to 43km in the map (equivalent to 21km on
> the
> ground) may arise."
>
> Therefore you would need to use the WKT EXTENSION PROJ4:
>
> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
> +k=1.0 +units=m +wktext +no_defs
>
> without +nadgrids=@null
>
> In GRASS, you could use this proj4 definition to create a new location,
> then
> import the data (the -o flag might be needed), then reproject to the
> desired
> CRS.
>
> Considering this particular CRS (EPSG:3857), it is strange than gdalwarp
> works while GRASS reprojection does not work because GRASS uses GDAL to
> convert WKT to proj4, then to GRASS terminology. Some debugging is
> needed to
> figure out why within GRASS, the conversion from WKT to proj4 (using
> GDAL)
> produces wrong results.

Thank you, yes that works. I looked at lib/proj/convert.c where I
assume the problem is. If I interpret it correctly it basically
discards +a and +b from the EXTENSION and instead picks WGS84
ellipsoid because of wkt

DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]]

But I don't really know what that extension actually means or how
GRASS should treat it.

This is the WKT which is processed:

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0],EXTENSION["PROJ4","+proj=merc
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +wktext +no_defs"]]

This WKT is processed correctly, GDAL converts this with OSRExportToProj4()
to

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +nadgrids=@null +wktext +no_defs

the only problem is +nadgrids=@null.

Unfortunately, GPJ_osr_to_grass() discards any datum and/or ellipsoid info
in the proj4 term at
https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/convert.c#L477

Datum and/or ellipsoid definitions are determined later on from the original
OGR spatial reference, which causes problems for EPSG:3857 and possibly
other CRS's.

Are there actually more cases like this 3857 we need to fear?

The reason for this special handling is that GRASS has its own datum and
ellipsoid definitions in lib/gis, i.e. in datum.table, datumtransform.table,
ellipse.table, ellipse.table.solar.system.

It seems like a review of GRASS handling of GDAL CRS definitions is
needed...

That sounds rather complicated. Would some workaround be an option?
Maybe checking if the datum from GDAL matches the datum from the proj4
string?

Anna

Markus M

Anna

>
> Markus M
>
>>
>> Hi Anna!
>>
>> Thank you very much for confirming that! I am working with the NAIP
>> imagery as well. :slight_smile:
>>
>> I have found that their original Geotiff assets work perfectly.
>>
>> In fact, I was very happy to stumble on to the fact that the complete
>> NAIP
>> archive (~250 terabytes) is available as a bucket drive on Amazon Web
>> Services (AWS). So I setup GRASS instances to process the tiles, then
>> download the processed, reprojected images compressed as JP2s. I am
>> paying
>> for the bandwidth and compute time, but I think its worth it for my
>> purposes. I’ll be able to process and download the imagery I need in
>> about
>> 60 days compared to over 300 days without AWS!
>>
>>
>> Cheers,
>>
>> Jeshua Lacock
>> Founder/Engineer
>> <3DTOPO.com>
>> GlassPrinted.com
>>
>> _______________________________________________
>> grass-user mailing list
>> grass-user@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/grass-user

13

On Fri, Sep 1, 2017 at 5:16 AM, Anna Petrášová <kratochanna@gmail.com> wrote:

On Thu, Aug 31, 2017 at 3:47 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová <kratochanna@gmail.com>
wrote:

On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua@3dtopo.com>
wrote:

On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna@gmail.com>
wrote:

On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com>
wrote:

To hopefully help troubleshoot; I just reprojected one of the raster
tiles (from epsg: 3857), into the source location of one of the
lonlat
vectors (reverse projections from my OP), and the datasets are
offset by the
same distances. Since the dimensions of the raster are being changed
(by
r.proj), it leads me to think it must be a datum or coordinate
system
misalignment and not a projection issue.

I have the same problem, working with NAIP imagery. It is related to:
https://trac.osgeo.org/grass/ticket/2229

I have to remove nadgrids: @null from the PROJ_INFO file to be able
to
reproject into that location, but then it is shifted. gdalwarp helps.

EPSG:3857 is problematic because it
“Uses spherical development of ellipsoidal coordinates. Relative to WGS
84 /
World Mercator (CRS code 3395) errors of 0.7 percent in scale and
differences in northing of up to 43km in the map (equivalent to 21km on
the
ground) may arise.”

Therefore you would need to use the WKT EXTENSION PROJ4:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +wktext +no_defs

without +nadgrids=@null

In GRASS, you could use this proj4 definition to create a new location,
then
import the data (the -o flag might be needed), then reproject to the
desired
CRS.

Considering this particular CRS (EPSG:3857), it is strange than gdalwarp
works while GRASS reprojection does not work because GRASS uses GDAL to
convert WKT to proj4, then to GRASS terminology. Some debugging is
needed to
figure out why within GRASS, the conversion from WKT to proj4 (using
GDAL)
produces wrong results.

Thank you, yes that works. I looked at lib/proj/convert.c where I
assume the problem is. If I interpret it correctly it basically
discards +a and +b from the EXTENSION and instead picks WGS84
ellipsoid because of wkt

DATUM[“D_WGS_1984”,SPHEROID[“WGS_1984”,6378137.0,298.257223563]]

But I don’t really know what that extension actually means or how
GRASS should treat it.

This is the WKT which is processed:

PROJCS[“WGS_1984_Web_Mercator_Auxiliary_Sphere”,GEOGCS[“GCS_WGS_1984”,DATUM[“D_WGS_1984”,SPHEROID[“WGS_1984”,6378137.0,298.257223563]],PRIMEM[“Greenwich”,0.0],UNIT[“Degree”,0.017453292519943295]],PROJECTION[“Mercator_Auxiliary_Sphere”],PARAMETER[“False_Easting”,0.0],PARAMETER[“False_Northing”,0.0],PARAMETER[“Central_Meridian”,0.0],PARAMETER[“Standard_Parallel_1”,0.0],PARAMETER[“Auxiliary_Sphere_Type”,0.0],UNIT[“Meter”,1.0],EXTENSION[“PROJ4”,“+proj=merc
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +wktext +no_defs”]]

This WKT is processed correctly, GDAL converts this with OSRExportToProj4()
to

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +nadgrids=@null +wktext +no_defs

the only problem is +nadgrids=@null.

Unfortunately, GPJ_osr_to_grass() discards any datum and/or ellipsoid info
in the proj4 term at
https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/convert.c#L477

Datum and/or ellipsoid definitions are determined later on from the original
OGR spatial reference, which causes problems for EPSG:3857 and possibly
other CRS’s.

Are there actually more cases like this 3857 we need to fear?

The reason for this special handling is that GRASS has its own datum and
ellipsoid definitions in lib/gis, i.e. in datum.table, datumtransform.table,
ellipse.table, ellipse.table.solar.system.

It seems like a review of GRASS handling of GDAL CRS definitions is
needed…

That sounds rather complicated. Would some workaround be an option?
Maybe checking if the datum from GDAL matches the datum from the proj4
string?

The proj4 string is also coming from GDAL (OSRExportToProj4()). In this case, the WKT definition contains a datum (WGS84), while the proj4 string does not contain a datum, instead a sphere defined by a and b. A possible workaround could be to use any datum/ellipsoid definition from the proj4 string and not from the OGR spatial reference, but that would introduce another CRS translation, i.e. another potential source of errors. It might be easier and safer to handle the special case EPSG:3857 or maybe more general the WKT attribute EXTENSION if it exists and contains a proj4 definition.

Markus M

Anna

Markus M

Hi Anna!

Thank you very much for confirming that! I am working with the NAIP
imagery as well. :slight_smile:

I have found that their original Geotiff assets work perfectly.

In fact, I was very happy to stumble on to the fact that the complete
NAIP
archive (~250 terabytes) is available as a bucket drive on Amazon Web
Services (AWS). So I setup GRASS instances to process the tiles, then
download the processed, reprojected images compressed as JP2s. I am
paying
for the bandwidth and compute time, but I think its worth it for my
purposes. I’ll be able to process and download the imagery I need in
about
60 days compared to over 300 days without AWS!

Cheers,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

14

On Fri, Sep 1, 2017 at 9:41 AM, Markus Metz <markus.metz.giswork@gmail.com> wrote:

On Fri, Sep 1, 2017 at 5:16 AM, Anna Petrášová <kratochanna@gmail.com> wrote:

On Thu, Aug 31, 2017 at 3:47 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová <kratochanna@gmail.com>
wrote:

On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua@3dtopo.com>
wrote:

On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna@gmail.com>
wrote:

On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua@3dtopo.com>
wrote:

To hopefully help troubleshoot; I just reprojected one of the raster
tiles (from epsg: 3857), into the source location of one of the
lonlat
vectors (reverse projections from my OP), and the datasets are
offset by the
same distances. Since the dimensions of the raster are being changed
(by
r.proj), it leads me to think it must be a datum or coordinate
system
misalignment and not a projection issue.

I have the same problem, working with NAIP imagery. It is related to:
https://trac.osgeo.org/grass/ticket/2229

I have to remove nadgrids: @null from the PROJ_INFO file to be able
to
reproject into that location, but then it is shifted. gdalwarp helps.

EPSG:3857 is problematic because it
“Uses spherical development of ellipsoidal coordinates. Relative to WGS
84 /
World Mercator (CRS code 3395) errors of 0.7 percent in scale and
differences in northing of up to 43km in the map (equivalent to 21km on
the
ground) may arise.”

Therefore you would need to use the WKT EXTENSION PROJ4:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +wktext +no_defs

without +nadgrids=@null

In GRASS, you could use this proj4 definition to create a new location,
then
import the data (the -o flag might be needed), then reproject to the
desired
CRS.

Considering this particular CRS (EPSG:3857), it is strange than gdalwarp
works while GRASS reprojection does not work because GRASS uses GDAL to
convert WKT to proj4, then to GRASS terminology. Some debugging is
needed to
figure out why within GRASS, the conversion from WKT to proj4 (using
GDAL)
produces wrong results.

Thank you, yes that works. I looked at lib/proj/convert.c where I
assume the problem is. If I interpret it correctly it basically
discards +a and +b from the EXTENSION and instead picks WGS84
ellipsoid because of wkt

DATUM[“D_WGS_1984”,SPHEROID[“WGS_1984”,6378137.0,298.257223563]]

But I don’t really know what that extension actually means or how
GRASS should treat it.

This is the WKT which is processed:

PROJCS[“WGS_1984_Web_Mercator_Auxiliary_Sphere”,GEOGCS[“GCS_WGS_1984”,DATUM[“D_WGS_1984”,SPHEROID[“WGS_1984”,6378137.0,298.257223563]],PRIMEM[“Greenwich”,0.0],UNIT[“Degree”,0.017453292519943295]],PROJECTION[“Mercator_Auxiliary_Sphere”],PARAMETER[“False_Easting”,0.0],PARAMETER[“False_Northing”,0.0],PARAMETER[“Central_Meridian”,0.0],PARAMETER[“Standard_Parallel_1”,0.0],PARAMETER[“Auxiliary_Sphere_Type”,0.0],UNIT[“Meter”,1.0],EXTENSION[“PROJ4”,“+proj=merc
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
+units=m +nadgrids=@null +wktext +no_defs”]]

This WKT is processed correctly, GDAL converts this with OSRExportToProj4()
to

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
+k=1.0 +units=m +nadgrids=@null +wktext +no_defs

the only problem is +nadgrids=@null.

Unfortunately, GPJ_osr_to_grass() discards any datum and/or ellipsoid info
in the proj4 term at
https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/convert.c#L477

Datum and/or ellipsoid definitions are determined later on from the original
OGR spatial reference, which causes problems for EPSG:3857 and possibly
other CRS’s.

Are there actually more cases like this 3857 we need to fear?

According to OGR SRS handling, yes.

The reason for this special handling is that GRASS has its own datum and
ellipsoid definitions in lib/gis, i.e. in datum.table, datumtransform.table,
ellipse.table, ellipse.table.solar.system.

It seems like a review of GRASS handling of GDAL CRS definitions is
needed…

That sounds rather complicated. Would some workaround be an option?
Maybe checking if the datum from GDAL matches the datum from the proj4
string?

The proj4 string is also coming from GDAL (OSRExportToProj4()). In this case, the WKT definition contains a datum (WGS84), while the proj4 string does not contain a datum, instead a sphere defined by a and b. A possible workaround could be to use any datum/ellipsoid definition from the proj4 string and not from the OGR spatial reference, but that would introduce another CRS translation, i.e. another potential source of errors. It might be easier and safer to handle the special case EPSG:3857 or maybe more general the WKT attribute EXTENSION if it exists and contains a proj4 definition.

Fixed in trunk r71523 where any PROJ4 string provided in a WKT EXTENSION is used. OGR SRS handling might also provide PROJ4_GRIDS which should also be considered somehow.

Markus M

Markus M

Anna

Markus M

Hi Anna!

Thank you very much for confirming that! I am working with the NAIP
imagery as well. :slight_smile:

I have found that their original Geotiff assets work perfectly.

In fact, I was very happy to stumble on to the fact that the complete
NAIP
archive (~250 terabytes) is available as a bucket drive on Amazon Web
Services (AWS). So I setup GRASS instances to process the tiles, then
download the processed, reprojected images compressed as JP2s. I am
paying
for the bandwidth and compute time, but I think its worth it for my
purposes. I’ll be able to process and download the imagery I need in
about
60 days compared to over 300 days without AWS!

Cheers,

Jeshua Lacock
Founder/Engineer
<3DTOPO.com>
GlassPrinted.com

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

15

On Sun, Oct 1, 2017 at 4:49 PM, Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Fri, Sep 1, 2017 at 9:41 AM, Markus Metz <markus.metz.giswork@gmail.com>
wrote:

On Fri, Sep 1, 2017 at 5:16 AM, Anna Petrášová <kratochanna@gmail.com>
wrote:
>
> On Thu, Aug 31, 2017 at 3:47 PM, Markus Metz
> <markus.metz.giswork@gmail.com> wrote:
> >
> >
> > On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová
> > <kratochanna@gmail.com>
> > wrote:
> >>
> >> On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz
> >> <markus.metz.giswork@gmail.com> wrote:
> >> >
> >> >
> >> > On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua@3dtopo.com>
> >> > wrote:
> >> >>
> >> >>
> >> >> > On Aug 27, 2017, at 8:54 PM, Anna Petrášová
> >> >> > <kratochanna@gmail.com>
> >> >> > wrote:
> >> >> >
> >> >> > On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock
> >> >> > <jeshua@3dtopo.com>
> >> >> > wrote:
> >> >> >>
> >> >> >> To hopefully help troubleshoot; I just reprojected one of the
> >> >> >> raster
> >> >> >> tiles (from epsg: 3857), into the source location of one of the
> >> >> >> lonlat
> >> >> >> vectors (reverse projections from my OP), and the datasets are
> >> >> >> offset by the
> >> >> >> same distances. Since the dimensions of the raster are being
> >> >> >> changed
> >> >> >> (by
> >> >> >> r.proj), it leads me to think it must be a datum or coordinate
> >> >> >> system
> >> >> >> misalignment and not a projection issue.
> >> >> >
> >> >> > I have the same problem, working with NAIP imagery. It is
> >> >> > related to:
> >> >> > https://trac.osgeo.org/grass/ticket/2229
> >> >> >
> >> >> > I have to remove nadgrids: @null from the PROJ_INFO file to be
> >> >> > able
> >> >> > to
> >> >> > reproject into that location, but then it is shifted. gdalwarp
> >> >> > helps.
> >> >
> >> > EPSG:3857 is problematic because it
> >> > "Uses spherical development of ellipsoidal coordinates. Relative to
> >> > WGS
> >> > 84 /
> >> > World Mercator (CRS code 3395) errors of 0.7 percent in scale and
> >> > differences in northing of up to 43km in the map (equivalent to
> >> > 21km on
> >> > the
> >> > ground) may arise."
> >> >
> >> > Therefore you would need to use the WKT EXTENSION PROJ4:
> >> >
> >> > +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
> >> > +y_0=0
> >> > +k=1.0 +units=m +wktext +no_defs
> >> >
> >> > without +nadgrids=@null
> >> >
> >> > In GRASS, you could use this proj4 definition to create a new
> >> > location,
> >> > then
> >> > import the data (the -o flag might be needed), then reproject to
> >> > the
> >> > desired
> >> > CRS.
> >> >
> >> > Considering this particular CRS (EPSG:3857), it is strange than
> >> > gdalwarp
> >> > works while GRASS reprojection does not work because GRASS uses
> >> > GDAL to
> >> > convert WKT to proj4, then to GRASS terminology. Some debugging is
> >> > needed to
> >> > figure out why within GRASS, the conversion from WKT to proj4
> >> > (using
> >> > GDAL)
> >> > produces wrong results.
> >>
> >> Thank you, yes that works. I looked at lib/proj/convert.c where I
> >> assume the problem is. If I interpret it correctly it basically
> >> discards +a and +b from the EXTENSION and instead picks WGS84
> >> ellipsoid because of wkt
> >>
> >> DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]]
> >>
> >> But I don't really know what that extension actually means or how
> >> GRASS should treat it.
> >>
> >>
> >> This is the WKT which is processed:
> >>
> >>
> >>
> >>
> >> PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0],EXTENSION["PROJ4","+proj=merc
> >> +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0
> >> +units=m +nadgrids=@null +wktext +no_defs"]]
> >
> > This WKT is processed correctly, GDAL converts this with
> > OSRExportToProj4()
> > to
> >
> > +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
> > +y_0=0
> > +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
> >
> > the only problem is +nadgrids=@null.
> >
> > Unfortunately, GPJ_osr_to_grass() discards any datum and/or ellipsoid
> > info
> > in the proj4 term at
> >
> > https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/convert.c#L477
> >
> > Datum and/or ellipsoid definitions are determined later on from the
> > original
> > OGR spatial reference, which causes problems for EPSG:3857 and
> > possibly
> > other CRS's.
>
> Are there actually more cases like this 3857 we need to fear?

According to OGR SRS handling, yes.

> >
> > The reason for this special handling is that GRASS has its own datum
> > and
> > ellipsoid definitions in lib/gis, i.e. in datum.table,
> > datumtransform.table,
> > ellipse.table, ellipse.table.solar.system.
> >
> > It seems like a review of GRASS handling of GDAL CRS definitions is
> > needed...
> >
>
> That sounds rather complicated. Would some workaround be an option?
> Maybe checking if the datum from GDAL matches the datum from the proj4
> string?

The proj4 string is also coming from GDAL (OSRExportToProj4()). In this
case, the WKT definition contains a datum (WGS84), while the proj4 string
does not contain a datum, instead a sphere defined by a and b. A possible
workaround could be to use any datum/ellipsoid definition from the proj4
string and not from the OGR spatial reference, but that would introduce
another CRS translation, i.e. another potential source of errors. It might
be easier and safer to handle the special case EPSG:3857 or maybe more
general the WKT attribute EXTENSION if it exists and contains a proj4
definition.

Fixed in trunk r71523 where any PROJ4 string provided in a WKT EXTENSION is
used. OGR SRS handling might also provide PROJ4_GRIDS which should also be
considered somehow.

This is great, I just tested it with NAIP imagery
(https://prd-tnm.s3.amazonaws.com/StagedProducts/NAIP/nc_2014/Imagery/m_3507811_se_17_1_20140618_20141118.jp2),
it works without problems. Not sure what else to test. Thank you!

Anna

Markus M

Markus M
>
>
> >>
> >>
> >> Anna
> >>
> >> >
> >> > Markus M
> >> >
> >> >>
> >> >> Hi Anna!
> >> >>
> >> >> Thank you very much for confirming that! I am working with the
> >> >> NAIP
> >> >> imagery as well. :slight_smile:
> >> >>
> >> >> I have found that their original Geotiff assets work perfectly.
> >> >>
> >> >> In fact, I was very happy to stumble on to the fact that the
> >> >> complete
> >> >> NAIP
> >> >> archive (~250 terabytes) is available as a bucket drive on Amazon
> >> >> Web
> >> >> Services (AWS). So I setup GRASS instances to process the tiles,
> >> >> then
> >> >> download the processed, reprojected images compressed as JP2s. I
> >> >> am
> >> >> paying
> >> >> for the bandwidth and compute time, but I think its worth it for
> >> >> my
> >> >> purposes. I’ll be able to process and download the imagery I need
> >> >> in
> >> >> about
> >> >> 60 days compared to over 300 days without AWS!
> >> >>
> >> >>
> >> >> Cheers,
> >> >>
> >> >> Jeshua Lacock
> >> >> Founder/Engineer
> >> >> <3DTOPO.com>
> >> >> GlassPrinted.com
> >> >>
> >> >> _______________________________________________
> >> >> grass-user mailing list
> >> >> grass-user@lists.osgeo.org
> >> >> https://lists.osgeo.org/mailman/listinfo/grass-user
> >