On Thu, Aug 31, 2017 at 3:47 PM, Markus Metz <[email protected]> wrote: > > > On Thu, Aug 31, 2017 at 5:59 AM, Anna Petrášová <[email protected]> > wrote: >> >> On Wed, Aug 30, 2017 at 4:46 PM, Markus Metz >> <[email protected]> wrote: >> > >> > >> > On Mon, Aug 28, 2017 at 9:11 PM, Jeshua Lacock <[email protected]> >> > wrote: >> >> >> >> >> >> > On Aug 27, 2017, at 8:54 PM, Anna Petrášová <[email protected]> >> >> > wrote: >> >> > >> >> > On Tue, Aug 22, 2017 at 7:46 PM, Jeshua Lacock <[email protected]> >> >> > 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? Anna > 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 >> >> [email protected] >> >> https://lists.osgeo.org/mailman/listinfo/grass-user > _______________________________________________ grass-user mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/grass-user
