[R-sig-Geo] Projection question...

[Barry Rowlingson]

On Thu, Nov 11, 2010 at 10:20 PM, Jonathan Greenberg
<greenberg at ucdavis.edu> wrote:
> Folks:
>
> Robert has a version of this email concerning the implementation in
> raster, but I'm wondering if there's something going on with rgdal or
> proj?  Can anyone think of a reason why reprojecting a file of North
> America with the following projection info:
>
> "+proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107 +x_0=5632642
> +y_0=4612546 +ellps=WGS84"
>
> to
>
> "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
>
> would be losing all data north of latitude = 46.63769 (even if there
> is clearly data above this latitude)?  Here's the raster info of the
> file I'm trying to reproject:
>
>> unsynced
> class       : RasterLayer
> filename    :
> nrow        : 277
> ncol        : 349
> ncell       : 96673
> min value   : 0
> max value   : 10.1305
> projection  : +proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107
> +x_0=5632642 +y_0=4612546 +ellps=WGS84
> xmin        : -16231.49
> xmax        : 11313351
> ymin        : -16231.5
> ymax        : 8976020
> xres        : 32462.99
> yres        : 32463
>
> A reprojection works in the sense that all data south of latitude =
> 46.63769 is placed in the right location.  Thoughts?  If there isn't
> any obvious answer from this non-code post, I'll append some code
> based on a new R package I've been developing so the specific datasets
> can be examined more closely.

 What do you mean by the raster losing data north of 46.6379 degrees?
And what are you using to reproject? (Because a raster normally needs
transforming to a new grid in the new CRS, spTransform wont do it).

 I see no problems doing an spTransform on a set of points in the
range of your raster summary:

pts=data.frame(x=seq(-16231.49,11313351,len=30),y=seq(-16231.5,8976020,len=30),z=1:30)
coordinates(pts)=~x+y
proj4string(pts)="+proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107
+x_0=5632642 +y_0=4612546 +ellps=WGS84"
spTransform(pts,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
+towgs84=0,0,0"))

gives me 30 points with no NA's or anything unexpected...

Barry