[GRASS-user] Re-projected Data Position Mismatch

Markus Metz markus.metz.giswork at gmail.com
Thu Aug 31 12:47:15 PDT 2017


On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová <kratochanna at gmail.com>
wrote:
>
> On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz
> <markus.metz.giswork at gmail.com> wrote:
> >
> >
> > On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <jeshua at 3dtopo.com>
wrote:
> >>
> >>
> >> > On Aug 27, 2017, at 8:54 PM, Anna Petrášová <kratochanna at gmail.com>
> >> > wrote:
> >> >
> >> > On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <jeshua at 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.  :)
> >>
> >> 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 at lists.osgeo.org
> >> https://lists.osgeo.org/mailman/listinfo/grass-user
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20170831/461f78d2/attachment.html>


More information about the grass-user mailing list