On Sun, Oct 1, 2017 at 4:49 PM, Markus Metz <[email protected]> wrote: > > > On Fri, Sep 1, 2017 at 9:41 AM, Markus Metz <[email protected]> > wrote: >> >> >> >> On Fri, Sep 1, 2017 at 5:16 AM, Anna Petrášová <[email protected]> >> wrote: >> > >> > 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? > > 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 >> > >> >> [email protected] >> > >> >> https://lists.osgeo.org/mailman/listinfo/grass-user >> > > >> _______________________________________________ grass-user mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/grass-user
