[R-sig-Geo] proj4string read by readOGR doesn't seem to precisely specify source shapefile projection

Glenn Stauffer gestauffer at gmail.com
Wed Jun 14 19:48:33 CEST 2017


Thanks for the link. That does shed some light, but I'm not sure I entirely understand. The .prj file does not store the EPSG ID, but it looks like the shapefile properties does include that information. Furthermore, when I overwrite the .prj file as Shaun suggested, the properties of the associated shapefile now does show " WKID: 3071 Authority: EPSG " whereas it did not previously. But aside from that, I'm  not sure why readOGR didn't identify the datum information in the first place. When I add "+datum=NAD83" to the proj4string, ArcGIS seems to accept it as the same CRS, even though the original datum was NAD83 HARN (which I don't know how to specify in R). But I am pretty new to this.... For now, overwriting the .prj with the desired one seems to work for me (and it is easy to do in a script).
Glenn

-----Original Message-----
From: R-sig-Geo [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Edzer Pebesma
Sent: Wednesday, June 14, 2017 11:12 AM
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] proj4string read by readOGR doesn't seem to precisely specify source shapefile projection

Also relevant to this question is this discussion/issue:

https://github.com/edzer/sfr/issues/380

On 14/06/17 17:14, Shaun Walbridge wrote:
> Glen,
> 
> Round-tripping projection information can be tricky, in part because 
> the normalization routines vary between different stacks, and what 
> counts as the same form is also implementation dependent. In this case 
> (round-tripping with ArcGIS), you may be better off using the 
> arcgisbinding package [1], which can do the conversion to and from 
> ArcGIS native data types (including some that aren't available via 
> OGR/GDAL). It's easy enough to install [2], and also lets you build 
> Geoprocessing tools which speak directly to R. I imagine there's also 
> a pure R solution to the normalization, but it would likely have side 
> effects for existing packages, and may not be easily solved.
> 
> If you know that your data is fine and want a simple solution, you 
> could also just overwrite the new shapefile’s .prj file with a copy of 
> the original.
> 
> Cheers,
> Shaun
> 
> 1. 
> https://urldefense.proofpoint.com/v2/url?u=https-3A__r-2Darcgis.github
> .io_assets_arcgisbinding.pdf&d=DwIGaQ&c=n6-cguzQvX_tUIrZOS_4Og&r=fCPRb
> 7QX-vd5bnO9gIJHCiX852SVUtyYX--xtCKtpfk&m=9su2o7yjMBqHmqZT3kZAlTRZwPC_W
> 0Tm65yvo7NmjLw&s=pIBxR3icnDLcn639TzBc0mkXy6Tv7G9Ip7xcWSAMweU&e=
> 2. 
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_R-2DAr
> cGIS_r-2Dbridge-2Dinstall&d=DwIGaQ&c=n6-cguzQvX_tUIrZOS_4Og&r=fCPRb7QX
> -vd5bnO9gIJHCiX852SVUtyYX--xtCKtpfk&m=9su2o7yjMBqHmqZT3kZAlTRZwPC_W0Tm
> 65yvo7NmjLw&s=JADt055JJUuCG1eltd7N9o5f_fu_geqGlnzO62SGtGM&e=
> 
> On 6/14/17, 9:35 AM, "R-sig-Geo on behalf of Glenn Stauffer" <r-sig-geo-bounces at r-project.org on behalf of gestauffer at gmail.com> wrote:
> 
>     I have a shapefile that I imported into R using readOGR from the rgdal
>     package. I do a little bit of work with it, like adding attribute
>     information, etc, then export it as an ESRI shapefile again, with a new
>     name. However, when I bring both the original and new shapefile into ArcGIS,
>     it tells me that the CRS do not match. 
>     
>     So, noting that all the projection parameters remain the same, but the
>     projection and coordinate system names are different, and the datum name is
>     dropped, my questions are:
>     
>     1.	Is the second CRS the same as the first?
>     2.	If so, why did the names change, and why does ArcGIS no longer
>     recognize it as the same?
>     3.	If not, how did it get changed?
>     4.	Can the proj4string be modified to be more specific, and if so, why
>     did readOGR not already do this to preserve all the information?
>     
>     I can use the new shapefiles just fine and readOGR interprets them as
>     identical, but it is a bit vexing that ArcGIS does not. And, I could of
>     course define it again in ArcGIS, but part of the motivation to work in R is
>     to obviate the need to point and click for many files.
>     
>     I appreciate any insights or enlightenment. 
>     
>     Thanks,
>     
>     Glenn
>     
>     Here is the original projection information from ArcGIS:
>     
>     Projected Coordinate System:    NAD_1983_HARN_Transverse_Mercator
>     
>     Projection: Transverse_Mercator
>     
>     False_Easting:  520000.00000000
>     
>     False_Northing: -4480000.00000000
>     
>     Central_Meridian:   -90.00000000
>     
>     Scale_Factor:   0.99960000
>     
>     Latitude_Of_Origin: 0.00000000
>     
>     Linear Unit:    Meter
>     
>      
>     
>     Geographic Coordinate System:   GCS_North_American_1983_HARN
>     
>     Datum:  D_North_American_1983_HARN
>     
>     Prime Meridian:     Greenwich
>     
>     Angular Unit:   Degree
>     
>     Here is the proj4string from R, which also agrees with the proj4string given
>     for this projection at www.spatialreference.org for epsg:3071 and also for
>     SR-ORG:7396.
>     
>     +proj=tmerc +lat_0=0 +lon_0=-90 +k=0.9996 +x_0=520000 +y_0=-4480000
>     +ellps=GRS80 +units=m +no_defs
>     
>     When I use writeOGR to export the SpatialPolygonsDataFrame with the above
>     proj4string, then bring it back into ArcGIS, the projection information is
>     given as the following, and is no longer recognized as identical to the
>     original.
>     
>     Projected Coordinate System:    Transverse_Mercator
>     
>     Projection: Transverse_Mercator
>     
>     false_easting:  520000.00000000
>     
>     false_northing: -4480000.00000000
>     
>     central_meridian:   -90.00000000
>     
>     scale_factor:   0.99960000
>     
>     latitude_of_origin: 0.00000000
>     
>     Linear Unit:    Meter
>     
>      
>     
>     Geographic Coordinate System:   GCS_GRS 1980(IUGG, 1980)
>     
>     Datum:  D_unknown
>     
>     Prime Meridian:     Greenwich
>     
>     Angular Unit:   Degree
>     
>      
>     
>      
>     
>     
>     	[[alternative HTML version deleted]]
>     
>     _______________________________________________
>     R-sig-Geo mailing list
>     R-sig-Geo at r-project.org
>     
> https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mail
> man_listinfo_r-2Dsig-2Dgeo&d=DwICAg&c=n6-cguzQvX_tUIrZOS_4Og&r=YFaRLkc
> UCdDkLrpTbNOUV9J1CwYBCTMwgm5tdQkRSm4&m=Ws5QzgjxjD1HGvdDia-3XoimuEiDzcY
> 4Wh-CLS905w0&s=gBdaQAI6IxVnnl8qdt_GTEeiM7f2_5dd15qzQOsQ_64&e=
>     
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

--
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/



More information about the R-sig-Geo mailing list