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

[Robert J. Hijmans]

Jonathan,

Hard to say what is going on. This example works for me:

r <- raster(nrow=277, ncol=349, xmn=-16231.49, xmx=11313351,
ymn=-16231.5, ymx=8976020, crs="+proj=lcc +lat_1=50 +lat_2=50
+lat_0=50 +lon_0=-107 +x_0=5632642 +y_0=4612546 +ellps=WGS84")

r[] = 1:ncell(r)

pr <- projectRaster(r, crs="+proj=longlat +ellps=WGS84 +datum=WGS84
+no_defs +towgs84=0,0,0", res=5)

Using a low resolution because this thing crossed the date line. In
any case, something like the below, where you specify the output
raster you want, is generally more sensible.

pr = raster(xmn=-168, xmx=-40, ymn=40, ymx=75)
res(pr) = 0.2
pr <- projectRaster(r, pr, progress='window')

Robert


On Fri, Nov 12, 2010 at 7:25 AM, Barry Rowlingson
<b.rowlingson at lancaster.ac.uk> wrote:
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>