<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
        {font-family:Calibri;
        panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
        {margin:0in;
        margin-bottom:.0001pt;
        font-size:11.0pt;
        font-family:"Calibri","sans-serif";}
a:link, span.MsoHyperlink
        {mso-style-priority:99;
        color:blue;
        text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
        {mso-style-priority:99;
        color:purple;
        text-decoration:underline;}
span.EmailStyle17
        {mso-style-type:personal-compose;
        font-family:"Calibri","sans-serif";
        color:windowtext;}
.MsoChpDefault
        {mso-style-type:export-only;
        font-family:"Calibri","sans-serif";}
@page WordSection1
        {size:8.5in 11.0in;
        margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
        {page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal">Hi,<o:p></o:p></p>
<p class="MsoNormal">I’ve encountered what appears to be a concurrency issue when calling OGRCoordinateTransformation::Transform from multiple threads.  What I’ve found is that when Transform is called from multiple asynchronous tasks, I’m getting inconsistent
 results given the same inputs.  This is pronounced with the projection used in the GINA Elevation dataset (Alaska); sample data can be found here:<o:p></o:p></p>
<p class="MsoNormal"><a href="http://ifsar.gina.alaska.edu/data/2013/DTM/IFSAR.SDMI.2013.FUGRO.DTM_N59W144/IFSAR.SDMI.2013.FUGRO.DTM_N59W144.tar.gz">http://ifsar.gina.alaska.edu/data/2013/DTM/IFSAR.SDMI.2013.FUGRO.DTM_N59W144/IFSAR.SDMI.2013.FUGRO.DTM_N59W144.tar.gz</a><o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">The projection used by this datasource is:<o:p></o:p></p>
<p class="MsoNormal">PROJCS["NAD_1983_CORS96_Alaska_Albers",GEOGCS["GCS_NAD_1983_CORS96",DATUM["NAD_1983_CORS96",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",55],PARAMETER["standard_parallel_2",65],PARAMETER["latitude_of_center",50],PARAMETER["longitude_of_center",-154],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Or in proj4:<o:p></o:p></p>
<p class="MsoNormal">+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon+0-154 +x_0=0 +y_0=0 +ellps=GRS80 + units=m +no_defs<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">The following code snippet reproduced the problem for me.  If I were comment out the async tasks and simply perform the transform synchronously, the outputs all remain consistent with the expected result.  Also make note of the CPLConfigOption
 “USE_PROJ_480_FEATURES”.  Setting this to “NO” will also produce consistent results, as indicated in the inline comments below.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">bool TestGdalConcurrency(const std::string& filename)<o:p></o:p></p>
<p class="MsoNormal">{<o:p></o:p></p>
<p class="MsoNormal">    auto dataset = GDALOpen(filename.c_str(), GDALAccess::GA_ReadOnly);<o:p></o:p></p>
<p class="MsoNormal">    auto projectionName = GDALGetProjectionRef(dataset);<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    // By default, this is set to yes<o:p></o:p></p>
<p class="MsoNormal">    // Configuring this setting to "NO" will force GDAL to acquire a lock when calling into proj<o:p></o:p></p>
<p class="MsoNormal">    // This is a workaround, as we'd prefer lockfree concurrent access to proj.<o:p></o:p></p>
<p class="MsoNormal">    CPLSetConfigOption("USE_PROJ_480_FEATURES", "YES");<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    // Setup "From" SpatialReference<o:p></o:p></p>
<p class="MsoNormal">    OGRSpatialReference wgs84SpatialReference;<o:p></o:p></p>
<p class="MsoNormal">    wgs84SpatialReference.SetWellKnownGeogCS("WGS84");<o:p></o:p></p>
<p class="MsoNormal">    // Setup "To" SpatialReference<o:p></o:p></p>
<p class="MsoNormal">    OGRSpatialReference spatialReference(projectionName);<o:p></o:p></p>
<p class="MsoNormal">    // Create transformation from WGS84 to the sourcefile's projection<o:p></o:p></p>
<p class="MsoNormal">    auto pFromPlanet = OGRCreateCoordinateTransformation(&wgs84SpatialReference, &spatialReference);<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    // Will also provide this consistent latitude, longitude as input<o:p></o:p></p>
<p class="MsoNormal">    // The expectation is that given constant inputs, we should always receive the exact same output.<o:p></o:p></p>
<p class="MsoNormal">    double longitude = -143.8996;<o:p></o:p></p>
<p class="MsoNormal">    double latitude = 60.0074;<o:p></o:p></p>
<p class="MsoNormal">    double x = longitude;<o:p></o:p></p>
<p class="MsoNormal">    double y = latitude;<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    // Our truth data (we're confident in this answer since there's no concurrent access into proj yet.<o:p></o:p></p>
<p class="MsoNormal">    pFromPlanet->Transform(1, &x, &y, nullptr);<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    // Asynchronously call Transform from as many threads as I can get, every result should match the truth data<o:p></o:p></p>
<p class="MsoNormal">    std::vector<std::future<bool>> results;<o:p></o:p></p>
<p class="MsoNormal">    for (int i = 0; i < 100000; i++)<o:p></o:p></p>
<p class="MsoNormal">    {<o:p></o:p></p>
<p class="MsoNormal">        results.emplace_back(std::async(std::launch::async, [longitude, latitude, pFromPlanet, x, y]()<o:p></o:p></p>
<p class="MsoNormal">        {<o:p></o:p></p>
<p class="MsoNormal">            double _x = longitude;<o:p></o:p></p>
<p class="MsoNormal">            double _y = latitude;<o:p></o:p></p>
<p class="MsoNormal">            <o:p></o:p></p>
<p class="MsoNormal">            pFromPlanet->Transform(1, &_x, &_y, nullptr);<o:p></o:p></p>
<p class="MsoNormal">            return _x == x && _y == y;<o:p></o:p></p>
<p class="MsoNormal">        }));<o:p></o:p></p>
<p class="MsoNormal">    }<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    // Sum the number of failed comparisons<o:p></o:p></p>
<p class="MsoNormal">    auto numFailures = 0;<o:p></o:p></p>
<p class="MsoNormal">    for (auto& result : results)<o:p></o:p></p>
<p class="MsoNormal">    {<o:p></o:p></p>
<p class="MsoNormal">        if (!result.get())<o:p></o:p></p>
<p class="MsoNormal">            ++numFailures;<o:p></o:p></p>
<p class="MsoNormal">    }<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    OGRCoordinateTransformation::DestroyCT(pFromPlanet);<o:p></o:p></p>
<p class="MsoNormal">    GDALClose(dataset);<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">    return numFailures == 0;<o:p></o:p></p>
<p class="MsoNormal">}<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">If there’s something I’m missing that would allow for lock-free coordinate transformations I’d be happy to hear about it.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thanks,<o:p></o:p></p>
<p class="MsoNormal">Alex<o:p></o:p></p>
</div>
</body>
</html>