shapefiles missing projection information

[Gavin Cotterill]

My apologies if some of my confusion isn't strictly ctmm-related or if this has been addressed elsewhere and I missed it.

Two weeks ago I had code humming along to create UD objects and write shapefiles. I got pulled into another project, and I assume what happened is that in the interim some packages were updated, (or system paths changed?) and now I'm having issues with rgdal, but really I do not know. My session info is tacked on at the end of this message. I was using R 3.6.3, but the issue persists after a fresh install today with R 4.0.4. I also tried reverting to an early 2020 rgdal version in 3.6.3, but that didn't help either.

In any case, the shapefiles are created, but whereas before my .prj files read as such:

PROJCS["WGS_1984_UTM_Zone_11N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

They now read:

GEOGCS["GCS_unknown",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Possibly this is related to proj4string -> proj6 changes? As discussed here:

https://gis.stackexchange.com/questions/374508/spproj4string-equivalent-in-the-proj6-framework

I used the same projection for all of my telemetry objects:

> projx <- as.character(sp::CRS("+init=epsg:4326 +proj=utm +zone=11")) 

> projection(`myTelemetryObject`) <- projx

And these carry over to my UD object:

> projection(UD2_08AS26) 

[1] "+init=epsg:4326 +proj=utm +zone=11 +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

But this issue exists whether I use writeShapefile on a UD object, eg:

# ctmm::writeShapefile(UD2_08AS26,

#                      "~/../Desktop/08AS26",

#                      file="AKDE_08AS26",

#                      level.UD=0.95,

#                      level=0.95)

or use rgdal::writeOGR on an `sp` SpatialPolygonsDataframe, eg:

# rgdal::writeOGR(testPolygons, 

#               dsn = "~/../Desktop/08AS26", 

#               layer = "AKDE_08AS26", 

#               driver = "ESRI Shapefile")

Only the .prj file appears different from my output two weeks ago. In fact, since I'm using the same projection for many individuals and had already written shapefiles, renaming a copy of an old .prj and putting in the 08AS26 directory gives me shapefiles that appear to work in Arc perfectly well. In other words, the only problem seems to be the creation of the .prj.

I suppose my questions, lacking a reprex, are: has anyone encountered this? If so, what did they learn? Also, is it possible to define the projection of a UD object that has already been created? If so it might help me troubleshoot this problem. I'd be happy to send a script and RDS file of a UD object along if someone is willing to look at this.

Any advice would be greatly appreciated!

Gavin 

Session info #1:

R version 3.6.3 (2020-02-29)

Platform: x86_64-w64-mingw32/x64 (64-bit)

Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:

[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   

[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          

[5] LC_TIME=English_United States.1252    

attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

 [1] forcats_0.5.0   stringr_1.4.0   dplyr_1.0.2     purrr_0.3.3    

 [5] readr_1.3.1     tidyr_1.0.2     tibble_2.1.3    ggplot2_3.3.3  

 [9] tidyverse_1.3.0 ctmm_0.5.10    

loaded via a namespace (and not attached):

 [1] Rcpp_1.0.6       cellranger_1.1.0 pillar_1.4.3     compiler_3.6.3  

 [5] dbplyr_2.0.0     tools_3.6.3      lubridate_1.7.4  jsonlite_1.7.2  

 [9] lifecycle_0.2.0  nlme_3.1-145     gtable_0.3.0     lattice_0.20-41 

[13] pkgconfig_2.0.3  rlang_0.4.10     reprex_0.3.0     cli_2.3.0       

[17] DBI_1.1.0        rstudioapi_0.11  haven_2.2.0      xfun_0.12       

[21] withr_2.4.1      xml2_1.2.2       httr_1.4.1       raster_3.4-5    

[25] fs_1.3.1         generics_0.0.2   vctrs_0.3.4      hms_0.5.3       

[29] grid_3.6.3       tidyselect_1.1.0 glue_1.4.2       R6_2.4.1        

[33] readxl_1.3.1     sp_1.4-5         modelr_0.1.6     magrittr_2.0.1  

[37] backports_1.1.5  scales_1.1.0     codetools_0.2-16 rvest_0.3.6     

[41] assertthat_0.2.1 colorspace_1.4-1 tinytex_0.20     stringi_1.4.6   

[45] munsell_0.5.0    broom_0.5.5      crayon_1.3.4    

Session info #2:

R version 4.0.4 (2021-02-15)

Platform: x86_64-w64-mingw32/x64 (64-bit)

Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:

[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   

[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          

[5] LC_TIME=English_United States.1252    

attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] ctmm_0.6.0

loaded via a namespace (and not attached):

 [1] compiler_4.0.4   rgdal_1.5-23     tools_4.0.4      rstudioapi_0.13 

 [5] sp_1.4-5         Rcpp_1.0.6       tinytex_0.29     raster_3.4-5    

 [9] codetools_0.2-18 grid_4.0.4       xfun_0.21        lattice_0.20-41 

[Christen Fleming]

Hi Gavin,

I'm interested in knowing what's going on here and if there's anything I need to change in the package for the proj6 update. I haven't changed anything in the package (regarding projections) recently, and I'm not aware that I need to. Nothing in sp is giving me any warnings.

One thing that's weird to me in your UD projection is "+towgs84=0,0,0". Doesn't that mean that something is being re-projected into WGS84?

When I run as.character(sp::CRS("+init=epsg:4326 +proj=utm +zone=11")), I don't get that extra part. I only get the datum and no_defs, which all makes sense to me.

When I look at UD objects that I make with my own projections, I also do not see anything like that.

If the above is not the issue... ctmm::writeShapefile.UD first runs ctmm::SpatialPolygonsDataFrame.UD and then rgdal::writeOGR. So I would next check if the sp object that you get out of ctmm::SpatialPolygonsDataFrame.UD is correct.

Best,

Chris

[Gavin Cotterill]

I think there's a few things going on, and surely user error is involved. I attached a reprex using CIlla.

Both of my current R installs ended up with the latest version of 'sp' somehow, and I'm reluctant to go back to an earlier version because there's enough driving me mad as-is. But I believe earlier versions of 'sp' accepted proj.4 arguments that later versions don't. The code that used to run for me without error:

> projx <- as.character(sp::CRS("+init=epsg:4326 +proj=utm +zone=11"))

> projection(myTelemObj) <- projx

now returns: 

Error in validate.projection(proj) : 

  A projected coordinate system must be specified.

That's straightforward enough, I just need to give it a valid proj.4 string.

> projx <- as.character(sp::CRS("+proj=utm +zone=11 +datum=WGS84"))

> projection(myTelemObj) <- projx

Runs without error. I'm not sure why it accepted the former previously. Which maybe it shouldn't have. EPSG 4326 isn't a projected coordinate system I think. The way these packages look up epsgs is just to sort through csvs in the library where the packages are installed, and those evidently change. So that could be part of it. I can illustrate this with sf, since I do have different versions of those installed:

Package version 'sf' 0.8.1, linking GEOS 3.6.1 GDAL 2.2.3 PROJ 4.9.3 does include epsg 42106 used in my reprex. Package version 'sf' 0.9.7, linking GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1 does not.

So (1) I was probably giving a bad/questionable proj4 string to start. (2) How packages reference coordinate systems changes over time. But how does that carry over to writeOGR? With both of my current installs, I can now redefine projections for telemetry objects that get carried over to my UD objects unchanged. But in writing them, my prj files end up with 'unknowns'. From my reprex using the two-point equidistant:

PROJCS["unknown",GEOGCS["GCS_unknown",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Two_Point_Equidistant"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Latitude_Of_1st_Point",-24.2850983757726],PARAMETER["Latitude_Of_2nd_Point",-25.024736901024],PARAMETER["Longitude_Of_1st_Point",31.7775407470691],PARAMETER["Longitude_Of_2nd_Point",31.8880063438124],UNIT["Meter",1.0]]

And when I redfined using epsg 42106:

PROJCS["unknown",GEOGCS["GCS_unknown",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",20.0],PARAMETER["Latitude_Of_Origin",5.0],UNIT["Meter",1.0]]

The polygons map directly on top of one another in Arc, so that's good. And the PROJCS argument at the beginning is present. But looks wrong?

Gavin

[Christen Fleming]

Yes, in ctmm you need to work in a projected coordinate system, in units of meters, and with minimal local distortion, which is what the default projection algorithm aims for. Otherwise, the units and calculations don't really make any sense. One day, far in the future, we might support ellipsoidal processes, but not any time soon.

I don't really know anything about prj files, but that string looks reasonable to me.

Best,

Chris