[GRASS-user] Re-projected Data Position Mismatch

[Anna Petrášová]

On Sun, Oct 1, 2017 at 4:49 PM, Markus Metz
<markus.metz.giswork at gmail.com> wrote:
>
>
> On Fri, Sep 1, 2017 at 9:41 AM, Markus Metz <markus.metz.giswork at gmail.com>
> wrote:
>>
>>
>>
>> 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?
>
> 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.  :)
>> > >> >>
>> > >> >> 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
>> > >
>>