1

I used the nsidcICE function located here: https://github.com/cran/raster/blob/master/R/nsidcICE.R and the following code to open a sea ice raster file. However, I would like to overlap vector points and "zoom in" on the raster to the Canadian Arctic (lat: 72 to 75, long: -90 to -80). I am new to R and GIS.

library(raster)
library(graticule)
library(rgdal)

tfile <- system.file("extdata",  "nt_20140728_f17_v1.1_n.bin", package = "graticule")
ice<-.rasterFromNSIDCFile("nt_20140728_f17_v1.1_n.bin")


meridians <- seq(-180, 160, by = 20)
parallels <- c(-80, -73.77, -68, -55, -45)
mlim <- c(-180, 180)
plim <- c(-88, -50)
grat <- graticule(lons = meridians, lats = parallels, xlim = mlim, ylim =     plim, proj = projection(ice))
labs <- graticule_labels(meridians, parallels, xline = -45, yline = -60, proj = projection(ice))
plot(ice, axes = FALSE)
plot(grat, add = TRUE, lty = 3)
text(labs, lab = parse(text= labs$lab), col= c("firebrick", "darkblue")[labs$islon + 1], cex = 0.85)
title(sprintf("Sea ice concentration %s", gsub(".bin", "", basename(ice))), cex.main = 0.8) 
title(sub = projection(ice), cex.sub = 0.6)
  • 1
    Welcome to GIS SE! As a new user, please be sure to take the tour to learn about this site's focused Q&A format. Then please edit your question to clarify where you are having trouble and on the results of any research you have done in an effort to try to answer the question. – Andy Dec 11 '18 at 14:44
  • You are brave using that function! I need to update it to allow use of later version files. BTW, you need parallels and plim that are positive. – mdsumner Dec 12 '18 at 19:52
0

try

plot(ice, axes = FALSE, xlim = c(-90,-80), ylim = c(72, 75))
plot(grat, add = TRUE, lty = 3, xlim = c(-90,-80), ylim = c(72, 75))

the xlim() and ylim() parameters allow you to customise the boundaries of your plot.

  • This can't work because the data is in polar stereographic, xlim and ylim don't map trivially to ranges in long/lat, and they are required units that are in metres (0, 0 is at the north pole for this file). – mdsumner Dec 12 '18 at 19:53
0

The NSIDC ice is in (north) Polar Stereographic, the easiest way to get a zoom is to use drawExtent() on the current plot, then use that to set up the plot

## <your plot above>
e <- drawExtent() ## interactive!
plot(e, asp = 1, xlab = "", ylab = "", axes = FALSE)
plot(ice, add = TRUE)
## the rest of <your plot above>

Note that I had to negate your parallels, they should be positive for the northern hemisphere (indicated by "_n" in the file name, and the "+lat_0" and "+lat_ts" in the projection(ice).

Alternatively, use projectExtent to discover the right xlim/ylim:

llex <- extent(-90, 80,  72, 75)
ppex <- projectExtent(raster(llex, crs = "+proj=longlat"), projection(ice))
plot(crop(ice), axes = FALSE)
## <the rest of your plot above>

I see

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.