Dear all,

Thanks for your help!
CRS() is working fine now:

 > library(rgdal)
 > My.Proj<-"+proj=tmerc +datum=WGS84 +lat_0=0 +lon_0=17 +x_0=600000 
+y_0=10000000"
 > CRS(My.Proj)
CRS arguments:
  +proj=tmerc +datum=WGS84 +lat_0=0 +lon_0=17 +x_0=600000 +y_0=10000000 
+ellps=WGS84 +towgs84=0,0,0

and the shapefiles I exported with writeORG() are properly projected. 
Roger, you are definitely right. The projection is +proj=tmerc, 
not+proj=utm. Thanks for the tips!
As for readOGR(), I have no clue about the problem but I got the 
information I needed by opening the *.prj file as Barry suggested so I 
don't need that function anymore.
Fortunately writeORG() is working well and saves the file in the right 
directory.

Thanks again!

Ivan

On 29.08.2013 13:26, Roger Bivand wrote:
> 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@r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>
>>
>


	[[alternative HTML version deleted]]

