[GRASS-user] Re-projected Data Position Mismatch

[Markus Metz]

On Fri, Sep 1, 2017 at 5:16 AM, Anna Petrášová <kratochanna at gmail.com>
wrote:
>
> On Thu, Aug 31, 2017 at 3:47 PM, Markus Metz
> <markus.metz.giswork at gmail.com> wrote:
> >
> >
> > 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.
>
>  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.  :)
> >> >>
> >> >> 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/20170901/cac7ab23/attachment.html>