0

I was assigned the task to clip a raster from .nc file from a .tif file.

edit (from comment): i want to extract temp. info from the .nc because i need to check the yearly mean temperature of a specific region. to be comparable the comparison has to occur on exactly the same area. The .nc file is larger than the previously checked area so i need to "clip" it to the extent of a .tif I have. The .tif data is in form 0|1 where it is 0 (or the .tif is smaller than the .nc) the .nc data should be "cliped". In the end i want to keep the .nc data but at the extent of the .tif while still retaining its resolution & projection. (.tif and .nc have different projections&pixel sizes)

Now ordinarily that wouldn't be a problem as i could use raster::crop. This doesn't deal with different projections and different pixel size/resolution though. (I still used it to generate an approximation, but it is not precise enough for the final infromation, as can be seen in the code snippet below). The obvious method to generate a more reliable dataset/rasterset would be to first use a method like raster::projectRaster or raster::sp.Transform # adding sp.transform was done in an edit to the original question and homogenize the datasets but this approach takes too much time, as i have to do this for quite a few .nc files.

I was told the best method would be to generate a normalized matrix from the smaller raster "clip_frame" and then just multiply it with the "nc_to_clip" raster. Doing so should prevent any errors through map projections or other factors. This makes a lot of sense to me in theory but I have no idea how to do this in practice. I would be very grateful to any kind of hint/code snippet or any other help.

I have looked at similar problems on StackOverflow (and other sites) like:

As I am not even sure how to frame the question correctly, I might have overlooked an answer to this problem, if so please point me there!

My (working) code so far, just to give you an idea of how I want to approach the topic (here using the crop-function).

#library(ncdf4)
library(raster)
library(rgdal)
library(tidyverse)

nc_list<-list.files(pattern = ".*0.nc$")                   # list of .nc files containing raster and temperature information
#nc_to_clip <- lapply(nc_list, raster, varname="GST")      # read in as raster
nc_to_clip < -raster(ABC.nc, vername="GST)

clip_frame <- raster("XYZ.tif")                         # read in .tif for further use as frame


mean_temp_from_raster<-function(input_clip_raster, input_clip_frame){  # input_clip_raster= raster to clip, input_clip_frame
  
  r2_coord<-rasterToPoints(input_clip_raster, spatial = TRUE)   # step 1 to extract coordinates
  map_clip <- crop(input_clip_raster, extent(input_clip_frame)) # use crop to cut the input_clip_raster (this being the function I have to extend on)
  temp<-raster::extract(map_clip, r2_coord@coords)              # step 2 to extract coordinates 
  temp_C<-temp*0.01-273.15                                      # convert kelvin*100 to celsius
  temp_C<-na.omit(temp_C)
  mean(temp_C)
  
  return_list<-list(map_clip, mean(temp_C))
  return(return_list)
  
}

mean_tempC<-lapply(nc_to_clip, mean_temp_from_raster,clip_frame)

Thanks!

PS: I don't have much experience working with .nc files and/or RasterLayers in R as I used to work with ArcGIS/Python (arcpy) for problems like this, which is not an option right now.

D.J
  • 1,180
  • 1
  • 8
  • 17
  • can you describe your goal more precisely? You say you need to `to clip a raster using another raster` but in your code you also mention resolution and you use `extract`, suggesting that your goal is different. – Robert Hijmans Aug 13 '20 at 09:24
  • i want to extract temp. info from the .nc because i need to check the yearly mean temperature of a specific region. to be comparable the comparison has to occur on exactly the same area. The .nc file is larger than the previously checked area so i need to "clip" it to the extent of a .tif I have. The .tif data is in form 0|1 where it is 0 (or the .tif is smaller than the .nc) the .nc data should be "cliped". In the end i want to keep the .nc data but at the extent of the .tif while still retaining its resolution & projection. (.tif and .nc have different projections&pixel sizes) – D.J Aug 13 '20 at 11:47
  • You could edit your question, rather than answering here. – Robert Hijmans Aug 13 '20 at 13:43

1 Answers1

0

Perhaps something like this?

library(raster)
nc <- raster(ABC.nc, vername="GST)
clip <- raster("XYZ.tif")

x <- as(extent(clip), "SpatialPolygons")
crs(x) <- crs(clip)
y <- sp::spTransform(x, crs(nc))

clipped <- crop(nc, y)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • hi, thanks for the answer! this works just fine but is not what i was looking for, as i didnt want to use either sp.transform or projectRaster (sorry for not specifying sp.transform as well as projectRaster). both take a lot of extra time to compute. it might be more akin to the raster calculator in qgis (with which i have no experience but if it is similar to arcgis might be a solution) i hoped for an "in-code" solution instead of relying on qgis for automatisation and as i use r doing stuff either in qgis or python is not necessarily ideal – D.J Aug 17 '20 at 06:33
  • In my anwer, `spTransform` takes about 0 seconds and you only need to do it once. I still do not understand your question – Robert Hijmans Aug 17 '20 at 08:51
  • first of all, your solution does, in outcome, answer my question. so im honestly grateful for the answer! my problem is, that i am supposed to do it without any projection transfromation. the idea is to take the clipping frame, keep the extent, set all coordinates to 0 or 1 depending on the value. then basically multiply it with the raster i want to "clip". this way i only have coordinates/values for the desired area. how to do this in practice is beyound me. the reason for it is precisely to avoid changing the projection. i have to transform hundreds of rasters so it does take quite long. – D.J Aug 17 '20 at 12:09
  • You *must* do some projection transformation at least one time, but once you have `clipped`, you can use it for all your other data? – Robert Hijmans Aug 17 '20 at 14:20
  • that is exactly the point. i was also taught to transform if i wanted to clip different projections. but... apparently there is a way to avoid transformation (in this case) which is exactly where i hit the wall. anyway - if i can't do it the way the boss wants, as long as it works reasonably close to expectations i'll work it out. the overall results don't vary too much between even the simple clip from my test-script and the spTransform solution so it shouldn't be of by a significant margin. – D.J Aug 18 '20 at 06:53