Re: [gdal-dev] OGRSpatialReference::IsSame returning TRUE unexpectedly
On mardi 5 septembre 2017 20:21:17 CEST Matthew Amato wrote: > Thanks for the detailed write up and quick fix! > > Just to clarify your last sentence, are you saying there are still > potential cases where IsSame returns TRUE when it shouldn't? Hopefully not. > I'm just wondering if I should ditch my optimization of checking > the spatial reference and always do a reproject (via VRT) even if it > results in a no-op in some cases. I assume the overhead of a no-op VRT is > fairly small and I'm concerned with correctness more than anything else. If "same for reprojection" is what matters for you, a safe way of checking is calling exportToProj4() on both SRS and checking for equality on the returned strings. -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Which Proj.4 transforms are available in GDAL?
Am 05.09.2017 um 20:50 schrieb Even Rouault: Apparently they also offer the datasets in WGS 84, so you can perhaps check which hypothesis is right from that. Both datasets align with the Custom CRS including 3-parms datum shift I included earlier in this thread. Greetings, André Joost ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGRSpatialReference::IsSame returning TRUE unexpectedly
Thanks for the detailed write up and quick fix! Just to clarify your last sentence, are you saying there are still potential cases where IsSame returns TRUE when it shouldn't? Or were you just emphasizing that this is a very tricky check to get right in all cases? I'm just wondering if I should ditch my optimization of checking the spatial reference and always do a reproject (via VRT) even if it results in a no-op in some cases. I assume the overhead of a no-op VRT is fairly small and I'm concerned with correctness more than anything else. Thanks again. On Tue, Sep 5, 2017 at 5:02 PM, Even Rouault wrote: > > I assume my understanding of IsSame is flawed so any clarification would > be > > > a big help. > > > > > > > Matt, > > > > Your understanding of IsSame() is correct. You just hit a bug. I've fixed > it per https://trac.osgeo.org/gdal/ticket/7029. You can see the changeset > code as a workaround to add in your own code if you need. > > > > IsSame() is a "smart" comparison, ie it tries to abstract away > unsignificant differences like the PROJCS text value that don't participate > really to the SRS and can vary over formats, data producers, etc... The > issue here is that web-mercator (EPSG:3857) and WGS84-mercator (EPSG:3395) > are represented the same way at the GEOGCS and projection levels. The only > difference is that GDAL creates a EXTENSION["PROJ4", ...] node in the > web-mercator case to really use the spherical definition. And IsSame() > before my fix neglected to compare that extension node. So you were a bit > unlucky since that must be pretty much the only pair of EPSG entries in > that situation. > > > > My initial tought was to compare the authority codes/nodes when they exist > on both sides, but I was wondering if there are not cases where you'd have > the same SRS registered in 2 different catalogs. Presumably within the same > catalog, a different code should mean a different CRS, but I'm not even > completely sure you couldn't find oddities. I was thinking to EPSG:3857 vs > the deprecated EPSG:3785, but this isn't a good example, since they differ > on the DATUM. So currently IsSame() returns FALSE on that pair, although I > guess one might consider them as equivalent (the deprecated version is just > a buggy expression on the current EPSG:3857). Or not... > > You also have the situation where really different CRS (IsSame() == FALSE) > ends up being translated the same way as proj.4 strings... > > > > So it is actually not obvious to define the semantics of IsSame(). It > should at least return FALSE for CRS that are clearly different (coordinate > transformation results in differences). > > > > Even > > > > -- > > Spatialys - Geospatial professional services > > http://www.spatialys.com > ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGRSpatialReference::IsSame returning TRUE unexpectedly
> I assume my understanding of IsSame is flawed so any clarification would be > a big help. > Matt, Your understanding of IsSame() is correct. You just hit a bug. I've fixed it per https:// trac.osgeo.org/gdal/ticket/7029. You can see the changeset code as a workaround to add in your own code if you need. IsSame() is a "smart" comparison, ie it tries to abstract away unsignificant differences like the PROJCS text value that don't participate really to the SRS and can vary over formats, data producers, etc... The issue here is that web-mercator (EPSG:3857) and WGS84-mercator (EPSG:3395) are represented the same way at the GEOGCS and projection levels. The only difference is that GDAL creates a EXTENSION["PROJ4", ...] node in the web-mercator case to really use the spherical definition. And IsSame() before my fix neglected to compare that extension node. So you were a bit unlucky since that must be pretty much the only pair of EPSG entries in that situation. My initial tought was to compare the authority codes/nodes when they exist on both sides, but I was wondering if there are not cases where you'd have the same SRS registered in 2 different catalogs. Presumably within the same catalog, a different code should mean a different CRS, but I'm not even completely sure you couldn't find oddities. I was thinking to EPSG:3857 vs the deprecated EPSG:3785, but this isn't a good example, since they differ on the DATUM. So currently IsSame() returns FALSE on that pair, although I guess one might consider them as equivalent (the deprecated version is just a buggy expression on the current EPSG:3857). Or not... You also have the situation where really different CRS (IsSame() == FALSE) ends up being translated the same way as proj.4 strings... So it is actually not obvious to define the semantics of IsSame(). It should at least return FALSE for CRS that are clearly different (coordinate transformation results in differences). Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] OGRSpatialReference::IsSame returning TRUE unexpectedly
Short version:
I don't understand why the below code snippet is returning TRUE. Unless
I'm missing something, EPSG 3857 is not the same as 3395, why is GDAL
reporting that they are?
auto epsg3857 = OGRSpatialReference();
epsg3857.importFromEPSG(3857);
auto epsg3395 = OGRSpatialReference();
epsg3395.importFromEPSG(3395);
if (epsg3857.IsSame(&epsg3395) == TRUE) {
std::cout << "Bug?" << std::endl;
}
Longer version:
I'm attempting to load in a GeoTIFF and detect whether it is already in
EPSG:3857 or needs to be reprojected. My function looks like this:
bool isSpatialReferenceSystem(GDALDataset &dataset, int epsg)
{
std::string in_srs_wkt = dataset.GetProjectionRef();
if (in_srs_wkt.length() == 0)
{
return false;
}
auto dataset_srs = OGRSpatialReference(in_srs_wkt.c_str());
auto target_srs = OGRSpatialReference();
target_srs.importFromEPSG(epsg);
return dataset_srs.IsSame(&target_srs) == TRUE;
}
This function worked as expected until a few days ago when I was given a
GeoTIFF with the following information (via gdalinfo):
Driver: GTiff/GeoTIFF
Files: data.tif
Size is 11779, 8832
Coordinate System is:
PROJCS["WGS 84 / World Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","3395"]]
Origin = (723576.68586748466,5474400.17280460610)
Pixel Size = (33.077567560570841,-33.077567560572312)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 723576.686, 5474400.173) ( 6d30' 0.00"E, 44d15' 0.34"N)
Lower Left ( 723576.686, 5182259.096) ( 6d30' 0.00"E, 42d20' 0.00"N)
Upper Right ( 1113197.354, 5474400.173) ( 10d 0' 0.08"E, 44d15' 0.34"N)
Lower Right ( 1113197.354, 5182259.096) ( 10d 0' 0.08"E, 42d20' 0.00"N)
Center ( 918387.020, 5328329.634) ( 8d15' 0.04"E, 43d17'57.55"N)
Band 1 Block=11779x1 Type=Byte, ColorInterp=Red
Mask Flags: PER_DATASET ALPHA
Band 2 Block=11779x1 Type=Byte, ColorInterp=Green
Mask Flags: PER_DATASET ALPHA
Band 3 Block=11779x1 Type=Byte, ColorInterp=Blue
Mask Flags: PER_DATASET ALPHA
Band 4 Block=11779x1 Type=Byte, ColorInterp=Alpha
Creating a GDALDataset from the image and calling my function returns TRUE
when I expect it to return false, after some debugging, I traced it back to
my simpler code snippet that shows 3857 and 3395 being treated as the same.
The odd thing is that if I tell GDAL to reproject the image to 3857 (via a
VRT), everything works. So GDAL clearly knows they are different under the
hood.
I assume my understanding of IsSame is flawed so any clarification would be
a big help.
Thanks,
Matt
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Which Proj.4 transforms are available in GDAL?
> Except that they have a scale factor with LCC 2SP, but it gets silently > dropped by GDAL/PROJ.4. Which is OK. The LCC 2SP formulation just needs the 2 standard parallels and the latitude of origin. Basically when GDAL imports a ESRI WKT and sees that there's 2 standard parallel, it considers it is a LCC_2SP. The scale factor is redundant there. proj.4 doesn't need it when the 2 standard parallels are provided; You need either the 2 standard parallels, or 1 standard parallel + the scale factor, and with some maths you can go from one form to the equivalent other one (the importFromESRI() method could probably be enhanced to check the consistency of the 4 parameters when present and emit a warning in case they are not consistent. Exercice left to contributors) Actually looking at proj.4 code, you must *not* include a +k_0 (other than 1) when using the +lat_1 +lat_2 formulation, otherwise the scale factor would be applied twice. > ... and that looks nothing similar to the Ocotepeque datum shift. They > are about 400 meters apart. Well that's expected if NAD27 and Octopeque are indeed 2 different datums: there must be a shift between. The EPSG database has also shifts from Octopeque to NAD27. Apparently they also offer the datasets in WGS 84, so you can perhaps check which hypothesis is right from that. -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Which Proj.4 transforms are available in GDAL?
Am 05.09.2017 um 19:43 schrieb Even Rouault: The shapefile clearly states DATUM["D_North_American_1927", and the download page (from the official agency) http://www.cnr.gob.sv/geoportal-cnr/ tells the same. Hum, actually looking at the full .prj, it looks similar to EPSG:5460 Yes, projection parameter values are the same, and the datum shift is not stored inside the Shapefile definition. See: $ gdalsrsinfo ESRI::dptoA_Lambert_NAD27.prj PROJ.4 : +proj=lcc +lat_1=13.316667 +lat_2=14.25 +lat_0=13.78 +lon_0=-89 +x_0=50 +y_0=295809.184 +datum=NAD27 +units=m +no_defs ... and that's the pitfall for Proj.4, resulting in zero datum shift because the grid shift files don't cover El Salvador. OGC WKT : PROJCS["IDGES_rev", GEOGCS["GCS_North_American_1927", DATUM["North_American_Datum_1927", SPHEROID["Clarke_1866",6378206.4,294.9786982]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Conformal_Conic_2SP"], PARAMETER["False_Easting",50.0], PARAMETER["False_Northing",295809.184], PARAMETER["Central_Meridian",-89.0], PARAMETER["Standard_Parallel_1",13.316667], PARAMETER["Standard_Parallel_2",14.25], PARAMETER["Scale_Factor",0.6704], PARAMETER["Latitude_Of_Origin",13.78], UNIT["Meter",1.0]] So beside the datum difference, they are equivalent. The .prj one uses a LCC2SP formulation of LCC, and the EPSG one a LCC1SP one. Except that they have a scale factor with LCC 2SP, but it gets silently dropped by GDAL/PROJ.4. So the NAD27 datum might be just an artifact of something in their production chain not being able to handle Octopeque 1935. Or just a matter of habit. My knowledge of El Salvador geodesy is rather limited ;-) ... or they keep a kind of mistery on it. El Salvador should have moved to SIRGAS-ES2007.8 by now, but I have not found any projected CRS based on that. Or if it is really NAD27 datum, I can indeed see that they are 3 NAD27->WGS84 transforms available for central America. $ ogrinfo pg:dbname=epsg -sql "select c.*, a.area_code, a.area_name from epsg_coordoperation c join epsg_area a on c.area_of_use_code = a.area_code where c.source_crs_code = 4267 and c.target_crs_code=4326 and area_south_bound_lat > 0 and area_north_bound_lat < 20 and area_west_bound_lon > -95 and area_east_bound_lon < -75" Two of them are specific of Panama. The remaining one is EPSG:1171 that applies to "Central America - Belize to Costa Rica" Its TOWGS84 params are 0.0, 125.0, 194.0, 0.0, 0.0, 0.0, 0.0 ... and that looks nothing similar to the Ocotepeque datum shift. They are about 400 meters apart. Greetings, André Joost ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Which Proj.4 transforms are available in GDAL?
> > The shapefile clearly states DATUM["D_North_American_1927", and the > download page (from the official agency) > http://www.cnr.gob.sv/geoportal-cnr/ tells the same. Hum, actually looking at the full .prj, it looks similar to EPSG:5460 See: $ gdalsrsinfo ESRI::dptoA_Lambert_NAD27.prj PROJ.4 : +proj=lcc +lat_1=13.316667 +lat_2=14.25 +lat_0=13.78 +lon_0=-89 +x_0=50 +y_0=295809.184 +datum=NAD27 +units=m +no_defs OGC WKT : PROJCS["IDGES_rev", GEOGCS["GCS_North_American_1927", DATUM["North_American_Datum_1927", SPHEROID["Clarke_1866",6378206.4,294.9786982]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Conformal_Conic_2SP"], PARAMETER["False_Easting",50.0], PARAMETER["False_Northing",295809.184], PARAMETER["Central_Meridian",-89.0], PARAMETER["Standard_Parallel_1",13.316667], PARAMETER["Standard_Parallel_2",14.25], PARAMETER["Scale_Factor",0.6704], PARAMETER["Latitude_Of_Origin",13.78], UNIT["Meter",1.0]] $ gdalsrsinfo EPSG:5460 PROJ.4 : +proj=lcc +lat_1=13.78 +lat_0=13.78 +lon_0=-89 +k_0=0.6704 +x_0=50 +y_0=295809.184 +ellps=clrk66 +towgs84=205,96,-98,0,0,0,0 +units=m +no_defs OGC WKT : PROJCS["Ocotepeque 1935 / El Salvador Lambert", GEOGCS["Ocotepeque 1935", DATUM["Ocotepeque_1935", SPHEROID["Clarke 1866",6378206.4,294.9786982138982, AUTHORITY["EPSG","7008"]], TOWGS84[205,96,-98,0,0,0,0], AUTHORITY["EPSG","1070"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","5451"]], PROJECTION["Lambert_Conformal_Conic_1SP"], PARAMETER["latitude_of_origin",13.78], PARAMETER["central_meridian",-89], PARAMETER["scale_factor",0.6704], PARAMETER["false_easting",50], PARAMETER["false_northing",295809.184], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["X",EAST], AXIS["Y",NORTH], AUTHORITY["EPSG","5460"]] So beside the datum difference, they are equivalent. The .prj one uses a LCC2SP formulation of LCC, and the EPSG one a LCC1SP one. One can easily check they are equivalent with: $ gdaltransform -s_srs "+proj=longlat +ellps=clrk66" -t_srs ESRI::dptoA_Lambert_NAD27.prj -88 14 608032.297052867 320003.359813003 0 $ gdaltransform -s_srs "+proj=longlat +ellps=clrk66" -t_srs EPSG:5460 -88 14 608032.297913122 320003.32045602 0 (the slight difference in northing might come from the .prj having a limited number of digits for Latitude_of_origin compared to the EPSG one) So the NAD27 datum might be just an artifact of something in their production chain not being able to handle Octopeque 1935. Or just a matter of habit. My knowledge of El Salvador geodesy is rather limited ;-) Or if it is really NAD27 datum, I can indeed see that they are 3 NAD27->WGS84 transforms available for central America. $ ogrinfo pg:dbname=epsg -sql "select c.*, a.area_code, a.area_name from epsg_coordoperation c join epsg_area a on c.area_of_use_code = a.area_code where c.source_crs_code = 4267 and c.target_crs_code=4326 and area_south_bound_lat > 0 and area_north_bound_lat < 20 and area_west_bound_lon > -95 and area_east_bound_lon < -75" Two of them are specific of Panama. The remaining one is EPSG:1171 that applies to "Central America - Belize to Costa Rica" Its TOWGS84 params are 0.0, 125.0, 194.0, 0.0, 0.0, 0.0, 0.0 Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Which Proj.4 transforms are available in GDAL?
Am 05.09.2017 um 16:56 schrieb Even Rouault: The problem is that there is no EPSG code for the projection, Ah, the -s_srs syntax is not limited to EPSG:. You can use a proj.4 string, inlined WKT, a filename that has WKT or proj.4 string in it, etc. See http://gdal.org/gdal_utilities.html and NAD27 datum is handled badly outside the US and Canada. The country in question is El Salvador, where they adopted NAD27, but it's not part of the NAD27 datum shift files. EPSG suggests a Ocotepeque datum, which is shifted to NAD27. Isn'it the datum used by ? $ gdalsrsinfo EPSG:5451 PROJ.4 : +proj=longlat +ellps=clrk66 +towgs84=205,96,-98,0,0,0,0 +no_defs OGC WKT : GEOGCS["Ocotepeque 1935", DATUM["Ocotepeque_1935", SPHEROID["Clarke 1866",6378206.4,294.9786982138982, AUTHORITY["EPSG","7008"]], TOWGS84[205,96,-98,0,0,0,0], AUTHORITY["EPSG","1070"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","5451"]] The TOWGS84 parameters are the ones of transformation http://epsg.io/6891 that is valid for El Savador (and neighbouring countries). There are other transformations from Octopeque 35 to WGS 84, but none have an area of use that include El Savador. They are also transformations from Octopeque 35 to NAD 27. The shapefile clearly states DATUM["D_North_American_1927", and the download page (from the official agency) http://www.cnr.gob.sv/geoportal-cnr/ tells the same. In fact, I have not found any Ocotepeque based GIS data available from El Salvador, only from the surrounding countries. Greetings, André Joost ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Which Proj.4 transforms are available in GDAL?
> The problem is that there is no EPSG code for the projection, Ah, the -s_srs syntax is not limited to EPSG:. You can use a proj.4 string, inlined WKT, a filename that has WKT or proj.4 string in it, etc. See http://gdal.org/gdal_utilities.html > and NAD27 > datum is handled badly outside the US and Canada. The country in > question is El Salvador, where they adopted NAD27, but it's not part of > the NAD27 datum shift files. EPSG suggests a Ocotepeque datum, which is > shifted to NAD27. Isn'it the datum used by ? $ gdalsrsinfo EPSG:5451 PROJ.4 : +proj=longlat +ellps=clrk66 +towgs84=205,96,-98,0,0,0,0 +no_defs OGC WKT : GEOGCS["Ocotepeque 1935", DATUM["Ocotepeque_1935", SPHEROID["Clarke 1866",6378206.4,294.9786982138982, AUTHORITY["EPSG","7008"]], TOWGS84[205,96,-98,0,0,0,0], AUTHORITY["EPSG","1070"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","5451"]] The TOWGS84 parameters are the ones of transformation http://epsg.io/6891 that is valid for El Savador (and neighbouring countries). There are other transformations from Octopeque 35 to WGS 84, but none have an area of use that include El Savador. They are also transformations from Octopeque 35 to NAD 27. Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Which Proj.4 transforms are available in GDAL?
Am 01.09.2017 um 21:05 schrieb Even Rouault: On vendredi 1 septembre 2017 20:42:34 CEST Andre Joost wrote: Am 01.09.2017 um 20:14 schrieb Even Rouault: ogr2ogr -s_srs EPSG:4326 -t_srs "+proj=lcc +lat_1=13.316667 +lat_2=14.25 +lat_0=13.78 +lon_0=-89 +x_0=50 +y_0=295809.184 +k_0=0.6704 +ellps=clrk66 +units=m +no_defs +towgs84=0,125,194,0,0,0,0 +wktext" LimitesNAD27.shp dptoA_WGS_1984.shp gdalsrsinfo LimitesNAD27.prj >>out.txt neither saves towgs84 or wktext Yes, a .prj file contains a ESRI flavour of WKT, which supports neither TOWGS84 or EXTENSION, so they get stripped. And more importantly the original EPSG code is also missing in ESRI WKT. So how is a datum shift outside of the EPSG code world supposed to be applied to a shapefile with ogr2ogr? You have to override the source SRS with -s_srs EPSG: if you know the code. One could potentially improve the situation by trying to match the .prj CRS name with the entries in the EPSG database to try recovering which code the ESRI WKT corresponds too (if there's a match) The problem is that there is no EPSG code for the projection, and NAD27 datum is handled badly outside the US and Canada. The country in question is El Salvador, where they adopted NAD27, but it's not part of the NAD27 datum shift files. EPSG suggests a Ocotepeque datum, which is shifted to NAD27. greetings, André Joost ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Return from Set/GetNoDataValue
Hi Andrew, > Perhaps this is a documentation issue, I think this is mostly an implementation issue in a niche driver. The changeset that introduced this ( https://trac.osgeo.org/gdal/changeset/21821 ) isn't very verbose of the motivation. I suspect that the intent was to benefit from PAM for a few things, but not wanting to generate a .aux.xml file when SetNoDataValue() is called. I've changed this in https://trac.osgeo.org/gdal/changeset/40003 to use the PAM .aux.xml file to laod non-default nodata value on GetNoDataValue(), and write into PAM on SetNoDataValue() if the value is different than the default one. Should be a reasonable compromise on the presumed original intent. Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Return from Set/GetNoDataValue
Perhaps this is a documentation issue, but I ran into:
double BTRasterBand::GetNoDataValue( int* pbSuccess /*= NULL */ )
{
if(pbSuccess != NULL)
*pbSuccess = TRUE;
return -32768;
}
CPLErr BTRasterBand::SetNoDataValue( double )
{
return CE_None;
}
It seems as if you can't set or fetch a nodata value, an error should be
returned. In this case, the code indicates no error. How is one supposed
to know (other than reading the code) that you can't set a value if success
is returned if this operation isn't supported?
Thanks,
--
Andrew Bell
andrew.bell...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Call for discussion on RFC68: C++11 compilation mode
Hi Tamas, this is great! Thanks Even > Hi Even, > > The new SDK-s have now been added. > > Best regards, > > Tamas > > 2017-08-29 15:18 GMT+02:00 Tamas Szekeres : > > Hi Even, > > > > Yes, I'm planning to include msvc2015/2017 versions shortly. > > > > Best regards, > > Tamas > > > > Sent from my iPhone > > > > 2017. aug. 14. dátummal, 19:02 időpontban Even Rouault < > > > > even.roua...@spatialys.com> írta: > > > Hi Tamas, > > > > > > Unless I'm mistaken your gisinternals.com builds don't include MSVC > > > > 2015 currently. Right ? Do you have any plans to do so ? Currently the > > GDAL > > AppVeyor target uses SDKSs from gisinternals to do builds with a number of > > third-party libraries. > > > > > Even > > > > > > > This is a call for discussion on RFC68 to set C++11 to be the minimum > > > > C++ > > > > > > language version for GDAL trunk. I've drafted the RFC based on work by > > > > Mateusz. > > > > > > > > https://trac.osgeo.org/gdal/wiki/rfc68_cplusplus11 > > > > > > > > There was initial discussions back in January of this year: > > > > > > > > https://lists.osgeo.org/pipermail/gdal-dev/2017-> > > > January/thread.html#45786 > > > > > > Cheers, > > > > -kurt > > > > Engineer at Google > > > > > > -- > > > Spatialys - Geospatial professional services > > > http://www.spatialys.com -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Tilesize in JPEG2000 and gdalinfo
On mardi 5 septembre 2017 09:20:14 CEST Rahkonen Jukka (MML) wrote: > This information may be less interesting once we have a released OpenJPEG > driver that can handle big tiles properly. Agreed. Sub-tile decoding with much more efficient RAM use is now implemented per https://github.com/uclouvain/openjpeg/pull/1010 . It is currently under review and should hopefully be merged into openjpeg master soon. A openjpeg 2.2.1 release with it should follow. I'll keep the list informed. In your use case, I see the image has RPCL progression order, which is, in theory, not so bad for efficient sub-window decoding when using network access (PCRL would probably be better here if doing only partial decoding at full resolution. RPCL allows for efficient sub-resolution decoding), but currently openjpeg ingests the whole codestream for the tiles it must decode totally or partially, so I/O wise, in a network context, this is not friendly for big single tiled images. I know some people would be interested in improving this. We'll see. With openjpeg + above mentionned pull request 1010, I tested gdal_translate /vsicurl/http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mara_v_25000_50/2013/N53/02m/1/N5144H.jp2 out.tif -srcwin 0 0 1024 1024 Took 1m42 in total, for just 4s of CPU time. So mostly the download time of 150 MB If you remove the /vsicurl/ prefix, then the image is completely downloaded in a single HTTP GET request (through the intermediate HTTP driver), which enables generally better throughput if you know that at the end the file will be entirely downloaded (so only for single tiled images). Note that this will consume 2 times the file size of RAM (once by the HTTP driver, and once by openjpeg) Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Tilesize in JPEG2000 and gdalinfo
Thanks, Even. This information may be less interesting once we have a released OpenJPEG driver that can handle big tiles properly. The reason for my question was this error message that I got from my colleague who is cutting excerpts from a big VRT Input file size is 1332094, 2304162 0...10...20...30...ERROR 1: Cannot decode tile, memory error ERROR 1: opj_decode() failed ERROR 1: /vsicurl/http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mara_v_25000_50/2013/N53/02m/1/N5144H.jp2, band 1: IReadBlock failed at X offset 0, Y offset 0: opj_decode() failed My guess is that error is due to this single tiled image. Most JPEG2000 images in the VRT are written with tilesize 1024x1024. -Jukka- Lähettäjä: Even Rouault [mailto:even.roua...@spatialys.com] Lähetetty: 5. syyskuuta 2017 12:09 Vastaanottaja: gdal-dev@lists.osgeo.org Kopio: Rahkonen Jukka (MML) Aihe: Re: [gdal-dev] Tilesize in JPEG2000 and gdalinfo Hi Jukka, > > If I run gdalinfo against the JPEG2000 image in > http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mara_v_25000_50/2013/N53/02m > /1/N5144H.jp2 the JP2ECW driver reports the blocksize as (256x256) and > JPEG2000OpenJPEG driver reports that the block size is (1024x1024). Yes for single tiled images, all JPEG2000 drivers will report a block size of much smaller dimension to avoid GDAL consuming too much memory. > > What interests me would be to know the fact that this image has only one > JPEG2000 tile with a tilesize of (11754x11754). Is there any > driver-independent way for getting this information? You can use the opj_dump utility from openjpeg, or a GDAL Python script https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/dump_jp2.py $ python swig/python/samples/dump_jp2.py /vsicurl/http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mara_v_25000_50/2013/N53/02m/1/N5144H.jp2 | grep "siz\"" 0 11754 <-- image width 11754 <-- image height 0 0 11754 <-- tile width 11754 <-- tile height 3 Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Tilesize in JPEG2000 and gdalinfo
Hi Jukka, > > If I run gdalinfo against the JPEG2000 image in > http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mara_v_25000_50/2013/N53/02m > /1/N5144H.jp2 the JP2ECW driver reports the blocksize as (256x256) and > JPEG2000OpenJPEG driver reports that the block size is (1024x1024). Yes for single tiled images, all JPEG2000 drivers will report a block size of much smaller dimension to avoid GDAL consuming too much memory. > > What interests me would be to know the fact that this image has only one > JPEG2000 tile with a tilesize of (11754x11754). Is there any > driver-independent way for getting this information? You can use the opj_dump utility from openjpeg, or a GDAL Python script https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/dump_jp2.py $ python swig/python/samples/dump_jp2.py /vsicurl/http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mara_v_25000_50/2013/N53/02m/1/N5144H.jp2 | grep "siz\"" 0 11754 <-- image width 11754 <-- image height 0 0 11754<-- tile width 11754<-- tile height 3 Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Tilesize in JPEG2000 and gdalinfo
Hi, If I run gdalinfo against the JPEG2000 image in http://kartat.kapsi.fi/files/orto/etrs-tm35fin/mara_v_25000_50/2013/N53/02m/1/N5144H.jp2 the JP2ECW driver reports the blocksize as (256x256) and JPEG2000OpenJPEG driver reports that the block size is (1024x1024). What interests me would be to know the fact that this image has only one JPEG2000 tile with a tilesize of (11754x11754). Is there any driver-independent way for getting this information? -Jukka Rahkonen- ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
