Projection question...
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