[R-sig-Geo] Converting State Plane Coordinates

[Jonathan Marc Bearak]

Hi,

Thanks for all the helpful comments I've received through this list.

Aside from solving the original issue, I've learned some incidentally useful things just by observing some of the code snippets.

As it turns out, I may have been calculating everything correctly -- but in a far less elegant way than by the methods suggested on this list, and I'm glad to know the proper way to handle this now regardless -- but making one fundamental error that hadn't occurred to me previously (as I'd thought it must have been my code, having no experience with GIS before now).  I had been checking, as a test case, to see whether a school was in a school attendance area using its geocoded address.  As I'm a Mac user, I quickly grabbed the test coordinates from geocoder.us, to avoid going back to the office, as all the GIS software seems to be Windows-based.  However, the coordinates that site generates are not sufficiently accurate, unlike Google Maps or Mappoint.

Regarding the datum, yes, it's NAD83.

Thanks,
Jonathan
Jonathan Marc Bearak
PhD Candidate in Sociology
Institute of Education Sciences Predoctoral Fellow
New York University
jonathan.bearak at nyu.edu




On 24 Aug, 2010,at 01:54 PM, Dan Putler <dan.putler at sauder.ubc.ca> wrote:

Roger,

This won't really matter since they are very close, but given the data 
source and a North America based company, my guess is the underlying 
datum is NAD83.

Dan

On 08/24/2010 10:45 AM, Roger Bivand wrote:
> On Tue, 24 Aug 2010, Alexandre Villers wrote:
>
> 
>> Hey,
>>
>> There is an ESRI code (ESRI:102318
>> <http://spatialreference.org/ref/esri/102318/>) corresponding to the
>> requested projection... However, I doubt CRS will take an ESRI code (Roger
>> ?).
>> 
> CRS("+init=esri:102318")
>
> but the ESRI version doesn't give a +datum=, so WGS84 will be assumed. The
> difference is the reversal of the +lat1= and +lat2= values, so probably
> not drastic.
>
> Roger
>
> 
>> Jonathan, have a look at the rgdal and sp packages help pages for the "How
>> To" in CRS()
>>
>> Best regards
>>
>> Alex
>>
>> Le 24/08/2010 17:58, Roger Bivand a écrit :
>> 
>>> On Tue, 24 Aug 2010, Alexandre Villers wrote:
>>>
>>> 
>>>> Hello,
>>>>
>>>> Have a look at spTransform (in rgdal package) and the EPSG code of the
>>>> desired projection (this is to me the easiest way not to mess with digits
>>>> and various copy errors that can be made while writing the projection
>>>> properties) at www.spatialreference.org
>>>>
>>>> then, you just have to do something like this
>>>>
>>>> library(rgdal)
>>>>
>>>> data<-read.table ("mydata.txt", h=T, sep=",") #your original dataset
>>>> coordinates(data)<-~ X + Y # where X and Y stand for the name of your
>>>> lat/lon columns
>>>> proj4string(data)<-CRS("+init=epsg:4326") #if your coordinates are in
>>>> WGS84
>>>> data.proj<-spTransfrom(data,
>>>> CRS("+init=epsg:the.correct.epsg.number.of.your.projection")
>>>> 
>>> Looks like:
>>>
>>> http://spatialreference.org/ref/epsg/32118/
>>>
>>> but has some different parameters. Search on "New York" in this site for
>>> alternatives.
>>>
>>> Roger
>>>
>>> 
>>>> HTH
>>>>
>>>> Alex
>>>>
>>>>
>>>> Le 24/08/2010 16:51, Jonathan Marc Bearak a écrit :
>>>> 
>>>>> Hi,
>>>>>
>>>>> I'm new to GIS and have been trying to convert latitude and longitude
>>>>> to/from state plane coordinates.
>>>>>
>>>>> I've tried using the project() program from the proj4 library to convert
>>>>> lat/lng to FIPS 3104 (New York State Long Island).
>>>>>
>>>>> No matter how I go about this, however, the coordinates come out wrong.
>>>>> E.g.,
>>>>> "+proj=lcc +lat_1=40.66666666666666
>>>>> +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000
>>>>> +y_0=0 +units=ft +no_defs +datum=NAD83",
>>>>> "+proj=lcc +a=6378137 +es=.0066943800229 +lon_0=-74
>>>>> +lat_1=41d2 +lat_2=40d40 +lat_0=40d10 +x_0=300000 +y_0=0 +units=ft
>>>>> +no_defs +datum=NAD83",
>>>>>
>>>>> E.g., if I try the inverse, to convert 1062791, 209606.6 to lat/lng,
>>>>> project() prints: 43.04762 -76.89626. The correct coordinates, however,
>>>>> are: 40.7416495 -073.7165681.
>>>>>
>>>>> I've been reading the mailing lists and searching Google and R project
>>>>> PDFs and manual pages without any luck for an embarrassingly long number
>>>>> of hours without any luck.
>>>>>
>>>>> Thanks for any help.
>>>>>
>>>>> Best,
>>>>> Jonathan
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at stat.math.ethz.ch
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>> 
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at stat.math.ethz.ch
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>> 
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>> 
>

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20100826/6a517ad0/attachment.html>