[GRASS-user] Re-projected Data Position Mismatch

[Anna Petrášová]

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"]]


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