title: PROJ4 to PROJ6 upgrade and "Discarded datum" warnings
body: <h1>Context</h1>
  <p>My questions are related to the changes induced by the upgrade from PROJ4 to PROJ6
  and the consequences in various R spatial packages (<code>sp</code>, <code>sf</code>, <code>raster</code>).</p>
  <p>We receive now a lot of warnings about “Discarded datum” that looks a bit worrying and I’m
  a little puzzled about what I should do with that. I can see that this can have dire consequences
  in some situations and that they can be ignored in other circumstances.</p>
  <p>It seems that I’m not the only one to be a little lost (see <a href="https://gis.stackexchange.com/questions/372692">here</a>).
  I hope that my questions with a specific reproducible example will help us to better understand the topic.</p>
  <p>I understand that we can <a href="https://gis.stackexchange.com/questions/371466">remove the warnings</a>,
  and I have read the context : <a href="https://www.r-spatial.org/r/2020/03/17/wkt.html" rel="noreferrer">r-spatial blog post</a>,
  <a href="https://cran.r-project.org/web/packages/rgdal/vignettes/PROJ6_GDAL3.html" rel="noreferrer">Migration to PROJ6/GDAL3</a>
  and <a href="https://rsbivand.github.io/ECS530_h19/ECS530_III.html" rel="noreferrer">these workshop notes</a> (mapview seems to
  handle this differently in more recent versions)</p>
  <h1>Questions</h1>
  <p><strong>Question 1 :</strong></p>
  <p>Probably a naive question :</p>
  <p>I understand that there is a need for a new notation/format (WKT) implemented
  in PROJ6 (eg because of a need of higher precision) but I don’t understand why there is
  a need to remove the datum part form the old proj4 string notation. Why not just
  keep it as it has always been (and implement the new features in the new WKT format/notation)</p>
  <p><strong>Question 2 :</strong></p>
  <p>It seems that we have 3 cases concerning the drop of the datum in the old proj4 format :</p>
  <ol>
  <li>no warning : the datum remains as in the old proj4string notation (<code>sf</code> défault ?)</li>
  <li>we have a warning “Discarded datum (…) but +towgs84= values preserved”</li>
  <li>we have a warning “Discarded datum (…)” (in that case the “+towgs84=” part of the
  string is discarded/dropped –&gt; <code>sp</code> default ?? )</li>
  </ol>
  <p>The example below illustrate different cases where we have these warnings.</p>
  <p>Why do we have these 3 different cases on the same CRS (here EPSG 31370)?<br />
  What are the consequences of the removal of the datum and/or <code>+towgs84</code> part ?<br />
  Should I be less worried by the second warning than by the third ?</p>
  <p><strong>Question 3 :</strong></p>
  <p>In the reproducible example below I try to extract the values from a raster corresponding
  to 3 points, the raster and the points having a different CRS.
  However depending on the approach used I get different results.
  I have the impression that this is related to this PROJ4 –&gt; PROJ6 changes and datum
  dropping but I might be wrong.
  I created this example only because I wanted to understand the consequences of these “datum dropping”
  in the crs</p>
  <p>I use the function <code>raster::extract</code> and 3 different general approaches (each time for both
  <code>sf</code> and <code>sp</code> objects for the points) from which I would expect the same output :</p>
  <ol>
  <li>no manual reprojection, I let <code>raster::extract</code> do the job to match the crs of the points and raster</li>
  <li>reprojecting manually the points toward the crs of the raster</li>
  <li>reprojecting manually the raster toward the crs of the points</li>
  </ol>
  <p>With the 3rd approach, I get 2 different set of values for <code>sp</code> or <code>sf</code> objects and
  I get yet a third set of values with the second (and first approach) (then the values are the same if
  I use <code>sp</code> or <code>sf</code> objects).</p>
  <ul>
  <li>Which values are the “correct” ones ? (probably <code>351.7868 236.4216 309.0073</code>)</li>
  <li>Why are some approaches failing (is this related to the warning messages / datum stuff ?) ?</li>
  <li>When I do a projection of a spatial object and I have these warnings, how can I check
  that my spatial objects have been “correctly” projected ?</li>
  </ul>
  <p><strong>Question 4 :</strong></p>
  <p>Is it possible to provide general recommendations about what should be done
  when we have these warning messages ?</p>
  <p>For example :</p>
  <ul>
  <li>If possible do not use the old proj4 string notation but use the WKT one instead
  (but this is not always possible. Here for example - if I’m not wrong I’m forced to use
  the old proj4 notation because this is what <code>raster</code> uses)</li>
  <li>If I knew the epsg code of the raster used here, I guess the safest appoach would be
  to reproject the points using that code with <code>sf</code> , eg : <code>st_transform(SF, crs = xxxx)</code></li>
  </ul>
  <hr />
  <h1>Reproducible example</h1>
  <h2>Create points and raster spatial data + check how the CRS is handled</h2>
  <h3>SF points spatial object</h3>
  <p>Briefly : the CRS is mainly stored in WKT format.
  The old proj4string is available on request and do not drop the datum/<code>towgs84</code> part</p>
  <pre class="lang-r prettyprint-override"><code>library(sf)
  #&gt; Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
  library(sp)
  library(raster)
  
  # create 3 points
  coo &lt;- data.frame(x = c(246000, 247000, 246500), y = c(184000, 186000, 185000))
  
  # create an sf spatial object
  SF &lt;- st_as_sf(coo, coords = c(&quot;x&quot;, &quot;y&quot;), crs = 31370)
  
  # Check the CRS : 
  # the proj4string includes the datum/+towgs84 information - no warning
  st_crs(SF)$proj4string
  #&gt; [1] &quot;+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs&quot;
  
  # Same value with `raster::crs` but with a 
  # Warning &quot;Discarded datum (...) but +towgs84= values preserved&quot;
  raster::crs(SF)
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded datum Unknown_based_on_International_1909_Hayford_ellipsoid in CRS definition,
  #&gt;  but +towgs84= values preserved
  #&gt; CRS arguments:
  #&gt;  +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
  #&gt; +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
  #&gt; +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs
  
  # WKT 
  st_crs(SF)
  #&gt; Coordinate Reference System:
  #&gt;   User input: EPSG:31370 
  #&gt;   wkt:
  #&gt; BOUNDCRS[
  #&gt;     SOURCECRS[
  #&gt;         PROJCRS[&quot;Belge 1972 / Belgian Lambert 72&quot;,
  #&gt;             BASEGEOGCRS[&quot;Belge 1972&quot;,
  #&gt;                 DATUM[&quot;Reseau National Belge 1972&quot;,
  #&gt;                     ELLIPSOID[&quot;International 1924&quot;,6378388,297,
  #&gt;                         LENGTHUNIT[&quot;metre&quot;,1]]],
  #&gt;                 PRIMEM[&quot;Greenwich&quot;,0,
  #&gt;                     ANGLEUNIT[&quot;degree&quot;,0.0174532925199433]],
  #&gt;                 ID[&quot;EPSG&quot;,4313]],
  #&gt;             CONVERSION[&quot;Belgian Lambert 72&quot;,
  #&gt;                 METHOD[&quot;Lambert Conic Conformal (2SP)&quot;,
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt;         ID[&quot;EPSG&quot;,1609],
  #&gt;         REMARK[&quot;Scale difference is given by information source as 0.999999. Given in this record in ppm to assist application usage. Very similar parameter values (to slightly less precision) used for BD72 to ETRS89: see code 1652.&quot;]]]
  cat(raster::wkt(SF)) # does not work with sf
  #&gt; Error in raster::wkt(SF): tentative d'obtenir le slot &quot;crs&quot; d'un objet (classe &quot;sf&quot;) qui n'est pas un objet S4
  </code></pre>
  <h3>SP points spatial object</h3>
  <p>Briefly : the CRS is mainly stored in proj4 string format and drops the datum and <code>towgs84</code> part
  unlike <code>sf</code>).
  the new WKT notation is stored as a comment in the CRS object but is different from <code>sf</code></p>
  <pre class="lang-r prettyprint-override"><code>SP &lt;- coo
  coordinates(SP) &lt;- ~x+y
  proj4string(SP) &lt;- CRS(&quot;+init=epsg:31370&quot;)
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded
  #&gt; datum Reseau_National_Belge_1972 in CRS definition
  
  # the proj4 string do not contain the `towgs84` part
  # Warning &quot;Discarded datum (...)&quot;
  CRS(&quot;+init=epsg:31370&quot;)
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded
  #&gt; datum Reseau_National_Belge_1972 in CRS definition
  #&gt; CRS arguments:
  #&gt;  +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
  #&gt; +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m
  #&gt; +no_defs
  # With `raster::crs` same proj4string but no warning
  raster::crs(SP)
  #&gt; CRS arguments:
  #&gt;  +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
  #&gt; +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m
  #&gt; +no_defs
  
  # WKT notation + a warning: this WKT is indeed different from the SF one (no datum here ?)
  cat(comment(CRS(&quot;+init=epsg:31370&quot;)))
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded
  #&gt; datum Reseau_National_Belge_1972 in CRS definition
  #&gt; PROJCRS[&quot;Belge 1972 / Belgian Lambert 72&quot;,
  #&gt;     BASEGEOGCRS[&quot;Belge 1972&quot;,
  #&gt;         DATUM[&quot;Reseau National Belge 1972&quot;,
  #&gt;             ELLIPSOID[&quot;International 1924&quot;,6378388,297,
  #&gt;                 LENGTHUNIT[&quot;metre&quot;,1]]],
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt;     USAGE[
  #&gt;         SCOPE[&quot;unknown&quot;],
  #&gt;         AREA[&quot;Belgium - onshore&quot;],
  #&gt;         BBOX[49.5,2.5,51.51,6.4]]]
  # Same output without warning
  cat(raster::wkt(SP))
  #&gt; PROJCRS[&quot;Belge 1972 / Belgian Lambert 72&quot;,
  #&gt;     BASEGEOGCRS[&quot;Belge 1972&quot;,
  #&gt;         DATUM[&quot;Reseau National Belge 1972&quot;,
  #&gt;             ELLIPSOID[&quot;International 1924&quot;,6378388,297,
  #&gt;                 LENGTHUNIT[&quot;metre&quot;,1]]],
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt;     USAGE[
  #&gt;         SCOPE[&quot;unknown&quot;],
  #&gt;         AREA[&quot;Belgium - onshore&quot;],
  #&gt;         BBOX[49.5,2.5,51.51,6.4]]]
  </code></pre>
  <h3>Raster</h3>
  <p>In <code>raster</code>, both the old proj4 notation and the new WKT notation seem to be stored.</p>
  <pre class="lang-r prettyprint-override"><code>r &lt;- raster::raster(system.file(&quot;external/test.grd&quot;, package=&quot;raster&quot;))
  raster::crs(r)
  #&gt; CRS arguments:
  #&gt;  +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889
  #&gt; +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs
  cat(raster::wkt(r))
  #&gt; PROJCRS[&quot;unknown&quot;,
  #&gt;     BASEGEOGCRS[&quot;unknown&quot;,
  #&gt;         DATUM[&quot;World Geodetic System 1984&quot;,
  #&gt;             ELLIPSOID[&quot;WGS 84&quot;,6378137,298.257223563,
  #&gt;                 LENGTHUNIT[&quot;metre&quot;,1]],
  #&gt;             ID[&quot;EPSG&quot;,6326]],
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt;         AXIS[&quot;(N)&quot;,north,
  #&gt;             ORDER[2],
  #&gt;             LENGTHUNIT[&quot;metre&quot;,1,
  #&gt;                 ID[&quot;EPSG&quot;,9001]]]]
  </code></pre>
  <p>This dataset do not contain <code>+towgs84</code> part in the proj4string. But when you do read
  a raster with a <code>+towgs84</code> part in the proj4string it seems to be dropped.<br />
  Non reproducible example :</p>
  <pre class="lang-r prettyprint-override"><code>GISfolder &lt;- &quot;/my/path&quot;
  tmp &lt;- raster(paste0(GISfolder, 'my_file.tif'))
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded
  #&gt; datum Unknown_based_on_International_1909_Hayford_ellipsoid in CRS definition
  raster::crs(tmp)
  #&gt; CRS arguments:
  #&gt;  +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=49.8333339
  #&gt; +lat_2=51.1666672333333 +x_0=150000.01256 +y_0=5400088.4378 +ellps=intl
  #&gt; +units=m +no_defs
  cat(raster::wkt(tmp))
  #&gt; PROJCRS[&quot;unknown&quot;,
  #&gt;     BASEGEOGCRS[&quot;unknown&quot;,
  #&gt;         DATUM[&quot;Unknown based on International 1909 (Hayford) ellipsoid&quot;,
  #&gt;             ELLIPSOID[&quot;International 1909 (Hayford)&quot;,6378388,297,
  #&gt;                 LENGTHUNIT[&quot;metre&quot;,1,
  #&gt;                     ID[&quot;EPSG&quot;,9001]]]],
  #&gt;         PRIMEM[&quot;Greenwich&quot;,0,
  #&gt;             ANGLEUNIT[&quot;degree&quot;,0.0174532925199433],
  #&gt;             ID[&quot;EPSG&quot;,8901]]],
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt;         AXIS[&quot;(N)&quot;,north,
  #&gt;             ORDER[2],
  #&gt;             LENGTHUNIT[&quot;metre&quot;,1,
  #&gt;                 ID[&quot;EPSG&quot;,9001]]]]
  </code></pre>
  <p>I should probably also explore what happens when we use the <code>stars</code> package instead of
  <code>raster</code> but this question
  is already quite long (and i have a lot of code build on the raster package)</p>
  <h2>Extract the values from the raster</h2>
  <h3>Using the automatic reprojection from <code>raster::extract</code></h3>
  <pre class="lang-r prettyprint-override"><code># extract the values from the raster, 
  # the function extract reprojects the points
  # in the same crs as the raster layer
  extract(r, SF)
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded datum Unknown_based_on_International_1909_Hayford_ellipsoid in CRS definition,
  #&gt;  but +towgs84= values preserved
  #&gt; Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
  #&gt; Raster
  #&gt; [1] 351.7868 236.4216 309.0073
  extract(r, SP)
  #&gt; Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
  #&gt; Raster
  #&gt; [1] 351.7868 236.4216 309.0073
  </code></pre>
  <h3>Reproject manually the points toward the crs of the raster</h3>
  SF with proj4string from the raster
  <pre class="lang-r prettyprint-override"><code>SF_proj &lt;- st_transform(SF, crs = raster::crs(r))
  extract(r, SF_proj)
  #&gt; [1] 351.7868 236.4216 309.0073
  </code></pre>
  SF with WKT from the raster
  <pre class="lang-r prettyprint-override"><code>SF_proj &lt;- st_transform(SF, crs = raster::wkt(r))
  extract(r, SF_proj)
  #&gt; [1] 351.7868 236.4216 309.0073
  </code></pre>
  SP with proj4string from the raster
  <pre class="lang-r prettyprint-override"><code>SP_proj &lt;- spTransform(SP, raster::crs(r))
  extract(r, SP_proj)
  #&gt; [1] 351.7868 236.4216 309.0073
  </code></pre>
  SP with WKT from the raster
  <p>The wkt fromat is not accepted by <code>sp::spTransform</code> –&gt; does not work</p>
  <pre class="lang-r prettyprint-override"><code># error
  SP_proj &lt;- spTransform(SP, raster::wkt(r))
  #&gt; Error in CRS(CRSobj): PROJ4 argument-value pairs must begin with +: PROJCRS[&quot;unknown&quot;,
  #&gt;     BASEGEOGCRS[&quot;unknown&quot;,
  #&gt;         DATUM[&quot;World Geodetic System 1984&quot;,
  #&gt;             ELLIPSOID[&quot;WGS 84&quot;,6378137,298.257223563,
  #&gt;                 LENGTHUNIT[&quot;metre&quot;,1]],
  #&gt;             ID[&quot;EPSG&quot;,6326]],
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt;         AXIS[&quot;(N)&quot;,north,
  #&gt;             ORDER[2],
  #&gt;             LENGTHUNIT[&quot;metre&quot;,1,
  #&gt;                 ID[&quot;EPSG&quot;,9001]]]]
  # extract(r, SP_proj)
  </code></pre>
  <h3>Reproject manually the raster toward the crs of the points</h3>
  Using the proj4string of SF (with datum)
  <p>–&gt; results different from the previous attempts</p>
  <pre class="lang-r prettyprint-override"><code># EPSG 31370 proj4 string with the datum:
  lambert72 &lt;- sf::st_crs(31370)$proj4string
  lambert72
  #&gt; [1] &quot;+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs&quot;
  # there is a warning when we project the raster but the full string seems to be used
  r2 &lt;- raster::projectRaster(r, crs = lambert72)
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded datum Unknown_based_on_International_1909_Hayford_ellipsoid in CRS definition,
  #&gt;  but +towgs84= values preserved
  raster::crs(r2)
  #&gt; CRS arguments:
  #&gt;  +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
  #&gt; +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
  #&gt; +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs
  
  extract(r2, SP)
  #&gt; [1] 341.6399 222.1028 301.2286
  </code></pre>
  Using the WKT of SF
  <p>–&gt; does not work because <code>raster::projectRaster</code> does not accept WKT format
  for its <code>crs</code> argument</p>
  <pre class="lang-r prettyprint-override"><code>lambert72 &lt;- sf::st_crs(31370)
  lambert72
  #&gt; Coordinate Reference System:
  #&gt;   User input: EPSG:31370 
  #&gt;   wkt:
  #&gt; BOUNDCRS[
  #&gt;     SOURCECRS[
  #&gt;         PROJCRS[&quot;Belge 1972 / Belgian Lambert 72&quot;,
  #&gt;             BASEGEOGCRS[&quot;Belge 1972&quot;,
  #&gt;                 DATUM[&quot;Reseau National Belge 1972&quot;,
  #&gt;                     ELLIPSOID[&quot;International 1924&quot;,6378388,297,
  #&gt;                         LENGTHUNIT[&quot;metre&quot;,1]]],
  #&gt;                 PRIMEM[&quot;Greenwich&quot;,0,
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt;             AREA[&quot;Belgium - onshore&quot;],
  #&gt;             BBOX[49.5,2.5,51.51,6.4]],
  #&gt;         ID[&quot;EPSG&quot;,1609],
  #&gt;         REMARK[&quot;Scale difference is given by information source as 0.999999. Given in this record in ppm to assist application usage. Very similar parameter values (to slightly less precision) used for BD72 to ETRS89: see code 1652.&quot;]]]
  r2 &lt;- raster::projectRaster(r, crs = lambert72)
  #&gt; Error in wkt(projto): tentative d'obtenir le slot &quot;crs&quot; d'un objet (classe &quot;crs&quot;) qui n'est pas un objet S4
  </code></pre>
  Using the proj4string of SP (without datum)
  <p>–&gt; results different from the previous attempts</p>
  <pre class="lang-r prettyprint-override"><code># EPSG 31370 proj4 string without the datum:
  lambert72 &lt;- sp::CRS(&quot;+init=epsg:31370&quot;)@projargs
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded
  #&gt; datum Reseau_National_Belge_1972 in CRS definition
  lambert72
  #&gt; [1] &quot;+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m +no_defs&quot;
  # warning
  r3 &lt;- raster::projectRaster(r, crs = lambert72)
  #&gt; Warning in showSRID(uprojargs, format = &quot;PROJ&quot;, multiline = &quot;NO&quot;): Discarded
  #&gt; datum Unknown_based_on_International_1909_Hayford_ellipsoid in CRS definition
  raster::crs(r3)
  #&gt; CRS arguments:
  #&gt;  +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
  #&gt; +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m
  #&gt; +no_defs
  
  extract(r3, coo)
  #&gt; [1] 348.5775 329.1199 277.2260
  </code></pre>
  <h1>sessionInfo</h1>
  <pre class="lang-r prettyprint-override"><code>sessionInfo()
  #&gt; R version 3.6.2 (2019-12-12)
  #&gt; Platform: x86_64-pc-linux-gnu (64-bit)
  #&gt; Running under: Ubuntu 18.04.4 LTS
  #&gt;
  #&gt; (...)
  #&gt; 
  #&gt; other attached packages:
  #&gt; [1] raster_3.3-13 sp_1.4-2      sf_0.9-5      knitr_1.29   
  #&gt; 
  #&gt; (...)
  #&gt; 
  </code></pre>
  <p>Used with :</p>
  <pre class="lang-r prettyprint-override"><code>GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
  </code></pre>
  <p><sup>Created on 2020-09-03 by the <a href="https://reprex.tidyverse.org" rel="noreferrer">reprex package</a> (v0.3.0)</sup></p>