[R-sig-Geo] [OT] Coordinates and pojection problems

[Roger Bivand]

On Thu, 19 Jun 2008, Paul Hiemstra wrote:

> Hi,
> You need to find out what the projection of the shapefile is, the 
> ogrinfo does not give this information for your file as you probably do 
> not have the .prj file. The R function spTransform needs the projection 
> info in a so called proj4string. The easiest way to define such a string 
> is through the epsg code, this code enables proj (the projection library 
> used in R) to find the parameter it needs. For example, the proj4string 
> for the Dutch national projection is:
>
> "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 
> +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m 
> +towgs84=565.04,49.91,465.84,-0.409394387,0.35970519561,-1.868491,4.0772 
> +no_defs"
>
> using epsg codes:
>
> "+init=epsg:28992"
>
> The bottomline is to find the projection of your second shapefile as 
> this is the most likely candidate to cause the difference. A good place 
> to ask questions regarding projections is the proj4-mailing list 
> (http://lists.maptools.org/mailman/listinfo/proj).

Right. In addition, do some digging into the history of coordinate 
reference systems for your country of interest, for example:

http://www.asprs.org/resources/grids/

has France on January 2001. I think that you may be in Nouvelle 
Triangulation Française, which has a Paris meridian:

EPSG <- make_EPSG()
EPSG[grep("NTF", EPSG$note), 1:2]

But that is guesswork. For enthusiasts, once you have a guess, use 
spTransform to geographical coordinates (CRS("+longlat +datum=WGS84")), 
and do writeOGR(obj, "obj.kml", "obj", driver="KML") of (a subset of) your 
object and see where it lands in Google Earth.

Roger

>
> hth and cheers,
> Paul
>
> Julien Barnier wrote:
>> Hi,
>>
>> I'm sorry to post a rather off-topic question here, but as a new sp
>> user and total novice in coordinates and projections questions, I'm
>> facing a problem you may help me to solve. Please ignore this post if
>> it is too much noise.
>>
>> I've got two map files (a MapInfo file and an ESRI shapefile) from two
>> different sources but which represent objects from the same french
>> region (one represent the cities, the other rivers and water
>> surfaces). 
>>
>> The problem is that when imported into R and plotted together, they
>> doesn't appear in the same place.  They seem to have the same
>> «longitude», but the «latitude» is different. If I plot them, the
>> first map ranges in the Y axis from 70,000 to 100,000 and the second
>> one from 2,060,000 to 2,140,000. The ranges along the x axis seem to
>> be the same, though.
>>
>> Here is the output of ogrinfo on the two files :
>>
>> ,----
>> | $ ogrinfo -so 69_iris.mid 69_iris
>> | Had to open data source read-only.
>> | INFO: Open of `69_iris.mid'
>> |       using driver `MapInfo File' successful.
>> | 
>> | Layer name: 69_iris
>> | Geometry: Unknown (any)
>> | Feature Count: 774
>> | Extent: (747910.000000, 2053380.000000) - (819740.000000,
>> | 2147622.300000)
>> | Layer SRS WKT:
>> | PROJCS[«unnamed»,
>> |     GEOGCS[«unnamed»,
>> |         DATUM[«MIF 9999,6,-168,-60,320,0,0,0,0,0»,
>> |             SPHEROID[«Clarke 1880»,6378249.145,293.465],
>> |             TOWGS84[-168,-60,320,-0,-0,-0,0]],
>> |         PRIMEM[«non-Greenwich»,0],
>> |         UNIT[«degree»,0.0174532925199433]],
>> |     PROJECTION[«Lambert_Conformal_Conic_2SP»],
>> |     PARAMETER[«standard_parallel_1»,45.90287723937],
>> |     PARAMETER[«standard_parallel_2»,47.69712276063],
>> |     PARAMETER[«latitude_of_origin»,46.8],
>> |     PARAMETER[«central_meridian»,2.337229104484],
>> |     PARAMETER[«false_easting»,600000],
>> |     PARAMETER[«false_northing»,2200000],
>> |     UNIT[«Meter»,1]]
>> | DepCom: String (5.0)
>> | Nom_Com: String (40.0)
>> | Iris: String (4.0)
>> | DComIris: String (9.0)
>> | Nom_Iris: String (43.0)
>> | Typ_Iris: String (1.0)
>> | Indic: String (1.0)
>> | Origine: String (1.0)
>> `----
>>
>> ,----
>> | $ ogrinfo -so Decoupage_Administratif/SURFDEAU.shp SURFDEAU
>> | INFO: Open of `Decoupage_Administratif/SURFDEAU.shp'
>> |       using driver `ESRI Shapefile' successful.
>> | 
>> | Layer name: SURFDEAU
>> | Geometry: Polygon
>> | Feature Count: 168
>> | Extent: (785619.190000, 68133.620000) - (818193.410000, 107443.340000)
>> | Layer SRS WKT:
>> | (unknown)
>> | NUMOBJ: String (11.0)
>> `----
>>
>> I've tried to play with some parameters in readShapePoly and
>> spTransform, but without success.
>>
>> Any help would be greatly appreciated.
>>
>> Thanks in advance,
>> 
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no