[R-sig-Geo] Import custom projection in R

[Roger Bivand]

On Thu, 29 Aug 2013, Pan wrote:

> Dear Roger,
>
> Thank you for your valuable suggestions. First, I tried to enter my.proj as 
> PROJ.4 string:
>
> > MyProj.PROJ.4<-"+proj=UTM +datum=WGS84 +lon_0=17 +x_0=600000 
> +y_0=10000000"
> >
> > library(rgdal)
> >
> > CRS("+proj=UTM +lon_0=17 +x_0=600000 +y_0=10000000")
> Error in CRS("+proj=UTM +lon_0=17 +x_0=600000 +y_0=10000000") :
>   unknown projection id
> >
> > CRS(paste("+proj=UTM +lon_0=17 +x_0=600000 +y_0=10000000"))
> Error in CRS(paste("+proj=UTM +lon_0=17 +x_0=600000 +y_0=10000000")) :
>   unknown projection id
>
> I tried to read a projected shapefile and retrieve the string from it because 
> I'm not particularly familiar with PROJ.4 and I was not sure
> whether the string was correct.

As Barry pointed out, +proj=utm is correct for Universal Transverse 
Mercator. But +proj=utm should be accompanied by +zone= and optionally 
+south in order to place the projection - you use +lon_0=, +x= and +y=, 
which I would expect to see rather used with +proj=tmerc, with +lat_0= to 
fix the northings.

> The function couldn't open the file:
>
> >
> > inputdir<-"C:/my work directory/"
> >
> > OGRSpatialRef(dsn=paste(inputdir), "my.proj")
> Error in OGRSpatialRef(dsn = paste(inputdir), "my.proj") :
>   Cannot open file
>
> although the directory is correct:
>
> > list.files(path=paste(inputdir), pattern="prj$")
> [1] "myproj.prj"

missing "."? No idea, maybe setwd("C:/my work directory/"), then dsn="."? 
As Barry said, open the offending file in notepad and scroll across to 
find the values.

Roger

> > list.files(path=paste(inputdir), pattern="^my.proj")
> [1] "my.proj.dbf" "my.proj.prj" "my.proj.sbn" "my.proj.sbx"
> "my.proj.shp" "my.proj.shx"
>
> The reason why readOGR() and similar functions can't read my file keeps being 
> a mystery but I will try to resolve it later.
> My current priority is to save the projection as object of the CRS-class and 
> have it available to transform coordinates, project spatial objects,
> etc. The on-line documentation led me always to the same point (i.e. enter it 
> as a PROJ.4 string) but I didn't manage to do it.
>
> Do you have any suggestion for me?
>
> Thanks in advance
>
> Ivan
>
> On 27.08.2013 09:53, Roger Bivand wrote:
>> On Tue, 27 Aug 2013, Pan wrote:
>> 
>>> Dear all,
>>> 
>>> I'm quite new to spatial data processing in R and I'm getting entangled in 
>>> projection systems, I would appreciate your support to find a solution. I 
>>> need to import my custom projection in R (my.proj), create projected 
>>> spatial objects (I usedSpatialPointsDataFrame {sp}) and export them as
>>> shapefiles (I usedwritePointsShape {maptools}). As I said, the first step 
>>> is particularly tricky to me!
>>> 
>>> I tried to import the projection both asmy.proj.shp andmy.proj.prj 
>>> usingreadOGR() {rgdal} but the result was just the same: it didn't work!
>>> 
>>> As you see below I've been trying different solutions to make the function 
>>> work but I still didn't understand what I'm supposed to assign to dsn and 
>>> layer:
>> 
>> No, reading (and re-reading) the help pages and their references may help. 
>> You have not explained why you need to read your "custom" projection. You 
>> could simply enter it as a PROJ.4 string. You appear to think that you 
>> should read a shapefile (including its optional *.prj file component) to 
>> retreive the string. If you don't need to do this, and can write it as a 
>> string in R, do that - CRS() of your string with rgdal loaded will check 
>> for obvious errors. If you have a representation that OGR recognises, such 
>> as a shapefile, you can say for example:
>> 
>> library(rgdal)
>> dsn <- system.file("vectors", package = "rgdal")[1]
>> # to access the vector examples provided
>> OGRSpatialRef(dsn, "cities")
>> 
>> Here, the dsn= argument is the directory containing the shapefile 
>> components, and cities is a layer; here, you can check that there are files 
>> ending in "prj", and for files starting with "cities":
>> 
>> list.files(path=dsn, pattern="prj$")
>> list.files(path=dsn, pattern="^cities")
>> 
>> Do check that the files are where you think they are, the most frequent 
>> cause of the error message you see being that they are somewhere else. Once 
>> you are done, use writeOGR(), not writePointShape() which is deprecated and 
>> which discards your coordinate reference system anyway, so defeating the 
>> point of your exercise.
>> 
>> Hope this clarifies,
>> 
>> Roger
>> 
>>> 
>>> inputdir<-"C:/my working directory/"
>>> 
>>> library(rgdal)
>>>
>>>  > Test<-readOGR(dsn=paste(inputdir),layer="my.proj")
>>> Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv 
>>> = use_iconv) :
>>>    Cannot open file
>>>
>>>  > Test<-readOGR(dsn=paste(inputdir,"my.proj",sep=""),layer="my.proj.shp")
>>> Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv 
>>> = use_iconv) :
>>>    Cannot open file
>>>
>>>  > Test<-readOGR(dsn="my.proj",layer="my.proj.shp")
>>> Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv 
>>> = use_iconv) :
>>>    Cannot open file
>>> 
>>> also the attempts withogrInfo()haven't been successful.
>>> 
>>> So I have two questions on this matter. How should I handle the issue 
>>> withreadOGR()? Is there any other way (possibly easy) to import a custom 
>>> projection in R?
>>> 
>>> Note: I'm using R 2.15.2 and rgdal 0.8-10.
>>> 
>>> Thanks in advance
>>> 
>>> Ivan
>>> 
>>>
>>>     [[alternative HTML version deleted]]
>>> 
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> 
>> 
>
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
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