Message-ID: <AANLkTikhMKmsUULqq2XXALr=wqWdyaJnz6dOnaZTEtTN@mail.gmail.com>
Date: 2010-11-13T06:39:52Z
From: Robert J. Hijmans
Subject: Projection question...
In-Reply-To: <AANLkTimaudR9jJ0=gxb6g=fZ4iQN2=wK3T-ygSCe4-6T@mail.gmail.com>
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
>