Как рассчитать% маскируемой области при применении функции маски?

Question

У меня есть два растровых объекта, обрезайте оба до одинаковой степени и замаскируйте все значения в растре 1, которые не соответствуют ie диапазону значений растра 2.

    suit_hab_45_50_Env <- futureEnv_45_50_cropped<maxValue(currentEnv_summer_masked)
    suit_hab_45_50_Env <- suit_hab_45_50_Env*futureEnv_45_50_cropped
    suit_hab_45_50_Env <- mask(futureEnv_45_50_cropped, suit_hab_45_50_Env, maskvalue=0)
    suit_hab_45_50_Env <- crop(suit_hab_45_50_Env, currentEnv_summer_masked)
    suit_hab_45_50_Env <- mask(suit_hab_45_50_Env, currentEnv_summer_masked)
    plot(suit_hab_45_50_Env)
    writeRaster(suit_hab_45_50_Env, filename = "suit_hab_45_50_Env", format = "GTiff", overwrite = TRUE)

Есть ли способ, которым R сообщает мне, сколько% площади растра 1 было замаскировано?

Или, другими словами: серый многоугольник = 100%, а перекрывающий растровый слой покрывает x% многоугольника.

Answer 1

Поскольку вы работаете в географической системе координат c, и особенно когда вы работаете в области высоких широт, вы не можете просто сравнить количество пикселей, отличных от NA, поскольку пиксели в более высоких широтах меньше, чем в более низкие широты. Вы должны использовать широту для вычисления площади каждого пикселя, а затем суммировать их, чтобы получить площадь каждого растра. Вот пример использования кадрированных и замаскированных растров Канады и функции raster::area, которая учитывает широту каждого пикселя при вычислении площади пикселей:

require(raster)
require(maptools)

##get shapefile of canada
data("wrld_simpl")
can_shp=wrld_simpl[which(wrld_simpl$NAME=="Canada"),]

##create empty raster of the world
rs=raster()

##crop canada
rs_can=crop(rs,can_shp)

##calaculate area of each pixel
crop_area=area(rs_can)
plot(crop_area) ##note cells are smaller at higher latitudes
lines(can_shp)

##calculate crop area
crop_area_total=sum(crop_area[])

##mask canada
mask_area=mask(crop_area,can_shp)
plot(mask_area)
lines(can_shp)

##calculate mask area
mask_area_total=sum(mask_area[],na.rm=T)
mask_area_total
# [1] 9793736 
# which is not far from the wikipedia reported 9,984,670 km2 
# the under-estimate is due to the coarse resolution of the raster

##calculate percentage
mask_area_total/crop_area_total*100
# [1] 48.67837

NB Вы неправильно пометили широту и долготу :)

Answer 2

Предлагаемый ответ работает нормально, однако сейчас я работаю с этим кодом. Работает отлично.

    suit_hab_45_50_temp_poly <- rasterToPolygons(suit_hab_45_50_temp, na.rm = TRUE)
    shapefile(suit_hab_45_50_temp_poly, filename="suit_hab_45_50_temp_poly.shp", 
    overwrite=TRUE)

    mosaic_summer_poly_arcmap <- spTransform(mosaic_summer_poly_arcmap, crs(suit_hab_45_50_temp_poly))

    mosaic_summer_poly_arcmap_area <- area(mosaic_summer_poly_arcmap)
    mosaic_summer_poly_arcmap_area_total <- sum(mosaic_summer_poly_arcmap_area[])
    mosaic_summer_poly_arcmap_area_total

    suit_hab_45_50_temp_area <- area(suit_hab_45_50_temp_poly)
    suit_hab_45_50_temp_area_total <- sum(suit_hab_45_50_temp_area[])
    suit_hab_45_50_temp_area_total

    hab_loss_45_50_temp_s <- 100- (suit_hab_45_50_temp_area_total/mosaic_summer_poly_arcmap_area_total*100)
    hab_loss_45_50_temp_s