<div dir="ltr"><div><div><div><div><div><div><div><br><br>On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová <<a href="mailto:kratochanna@gmail.com">kratochanna@gmail.com</a>> wrote:<br>><br>> On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz<br>> <<a href="mailto:markus.metz.giswork@gmail.com">markus.metz.giswork@gmail.com</a>> wrote:<br>> ><br>> ><br>> > On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <<a href="mailto:jeshua@3dtopo.com">jeshua@3dtopo.com</a>> wrote:<br>> >><br>> >><br>> >> > On Aug 27, 2017, at 8:54 PM, Anna Petrášová <<a href="mailto:kratochanna@gmail.com">kratochanna@gmail.com</a>><br>> >> > wrote:<br>> >> ><br>> >> > On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <<a href="mailto:jeshua@3dtopo.com">jeshua@3dtopo.com</a>><br>> >> > wrote:<br>> >> >><br>> >> >> To hopefully help troubleshoot; I just reprojected one of the raster<br>> >> >> tiles (from epsg: 3857), into the source location of one of the lonlat<br>> >> >> vectors (reverse projections from my OP), and the datasets are offset by the<br>> >> >> same distances. Since the dimensions of the raster are being changed (by<br>> >> >> r.proj), it leads me to think it must be a datum or coordinate system<br>> >> >> misalignment and not a projection issue.<br>> >> ><br>> >> > I have the same problem, working with NAIP imagery. It is related to:<br>> >> > <a href="https://trac.osgeo.org/grass/ticket/2229">https://trac.osgeo.org/grass/ticket/2229</a><br>> >> ><br>> >> > I have to remove nadgrids: @null from the PROJ_INFO file to be able to<br>> >> > reproject into that location, but then it is shifted. gdalwarp helps.<br>> ><br>> > EPSG:3857 is problematic because it<br>> > "Uses spherical development of ellipsoidal coordinates. Relative to WGS 84 /<br>> > World Mercator (CRS code 3395) errors of 0.7 percent in scale and<br>> > differences in northing of up to 43km in the map (equivalent to 21km on the<br>> > ground) may arise."<br>> ><br>> > Therefore you would need to use the WKT EXTENSION PROJ4:<br>> ><br>> > +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0<br>> > +k=1.0 +units=m +wktext  +no_defs<br>> ><br>> > without +nadgrids=@null<br>> ><br>> > In GRASS, you could use this proj4 definition to create a new location, then<br>> > import the data (the -o flag might be needed), then reproject to the desired<br>> > CRS.<br>> ><br>> > Considering this particular CRS (EPSG:3857), it is strange than gdalwarp<br>> > works while GRASS reprojection does not work because GRASS uses GDAL to<br>> > convert WKT to proj4, then to GRASS terminology. Some debugging is needed to<br>> > figure out why within GRASS, the conversion from WKT to proj4 (using GDAL)<br>> > produces wrong results.<br>><br>> Thank you, yes that works. I looked at lib/proj/convert.c where I<br>> assume the problem is. If I interpret it correctly it basically<br>> discards +a and +b from the EXTENSION and instead picks WGS84<br>> ellipsoid because of wkt<br>><br>> DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]]<br>><br>> But I don't really know what that extension actually means or how<br>> GRASS should treat it.<br>><br>><br>> This is the WKT which is processed:<br>><br>><br>> 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<br>> +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0<br>> +units=m +nadgrids=@null +wktext +no_defs"]]<br><br></div>This WKT is processed correctly, GDAL converts this with OSRExportToProj4() to<br><br>+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<br><br></div>the only problem is +nadgrids=@null.<br><br></div>Unfortunately, GPJ_osr_to_grass() discards any datum and/or ellipsoid info in the proj4 term at<br><a href="https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/convert.c#L477">https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/convert.c#L477</a><br><br></div>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.<br><br></div>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.<br><br></div>It seems like a review of GRASS handling of GDAL CRS definitions is needed...<br><br></div>Markus M<br><div><div><div><div><div><div><div>><br>><br>> Anna<br>><br>> ><br>> > Markus M<br>> ><br>> >><br>> >> Hi Anna!<br>> >><br>> >> Thank you very much for confirming that! I am working with the NAIP<br>> >> imagery as well.  :)<br>> >><br>> >> I have found that their original Geotiff assets work perfectly.<br>> >><br>> >> In fact, I was very happy to stumble on to the fact that the complete NAIP<br>> >> archive (~250 terabytes) is available as a bucket drive on Amazon Web<br>> >> Services (AWS). So I setup GRASS instances to process the tiles, then<br>> >> download the processed, reprojected images compressed as JP2s. I am paying<br>> >> for the bandwidth and compute time, but I think its worth it for my<br>> >> purposes. I’ll be able to process and download the imagery I need in about<br>> >> 60 days compared to over 300 days without AWS!<br>> >><br>> >><br>> >> Cheers,<br>> >><br>> >> Jeshua Lacock<br>> >> Founder/Engineer<br>> >> <3DTOPO.com><br>> >> GlassPrinted.com<br>> >><br>> >> _______________________________________________<br>> >> grass-user mailing list<br>> >> <a href="mailto:grass-user@lists.osgeo.org">grass-user@lists.osgeo.org</a><br>> >> <a href="https://lists.osgeo.org/mailman/listinfo/grass-user">https://lists.osgeo.org/mailman/listinfo/grass-user</a><br><br></div></div></div></div></div></div></div></div>