• 日期: 2027-07-20
  • 标签: 常规

Creating Transformation Matrices for Extracting Country-Level Climate Data

This script outlines the process of creating transformation matrices to extract country-level climate data from a global raster dataset. The script utilizes land cover shares or correlations as criteria for defining weights within countries.

1. Setup and Environment

# Clean up workspace
rm(list=ls())

# Clean and load packages
wants <- c('raster','rgdal','rworldmap','maptools','RColorBrewer','OIdata',
             'ncdf4','mFilter','plyr','lubridate','rworldxtra','Matrix','Cairo')
needs <- wants[!(wants %in% installed.packages()[,'Package'])]
if(length(needs)) install.packages(needs)
lapply(wants, function(i) require(i, character.only=TRUE))
rm(needs,wants)

# Directories
dir <- list()
dir$root <- dirname(getwd())
dir$data <- paste(dir$root,'/data',sep='')
dir$map <- paste(dir$root,'/data2/globalmap',sep='') 
dir$mask  <- paste(dir$root,'/data/Princeton/elevation_0.25deg.nc',sep='') 
dir$source <- paste(dir$root,'/data/Princeton/v3/0.25deg/daily',sep='')
dir$figures <- paste(dir$root,'/figures/SM',sep='') 
dir$weights <- paste(dir$root,'/data2/weights',sep='') 
lapply(dir, function(i) dir.create(i, recursive = T, showWarnings = F))

# Misc
misc <- list()
# Criteria for country-level extraction
# cropland: use cropland shares as weights within country
# pasture: use pasture shares as weights within country
# both: use cropland+pasture shares as weights within country
# cor: use significant correlations as indicator to include/exclude from weights
misc$criteria <- c('cropland','pasture','both') # 'cor'
#misc$criterianame <- list(cropland='Cropland', pasture='Pasture',both='Cropland + Pasture', cor='Correlation')
misc$criterianame <- list(cropland='Cropland', pasture='Pasture',both='Cropland + Pasture')
# Correlation criteria
misc$seasons <- c('12_1_2','3_4_5','6_7_8','9_10_11') # seasons to be considered
misc$varscrit <- c('tmax','tmin','prcp') # variables to be considered
misc$pval <- .05 # critical value of correlation to consider

2. Data Preparation and Preprocessing

# 1. Get raster data mask
mask <- r <- raster(dir$mask)
mask[] <- 1
mask[is.na(r[])] <- 0

# 2. Get polygon data
map <- readRDS(paste(dir$map,'/globalmap.RDS',sep=''))
map <- spTransform(map, CRS(projection(mask)))
map <- map[-which(map@data$NAME=='Antarctica'),] # dropAntarctica

# 3. Crop mask to exactly the climate data extend (the elevation file included the poles)
r <- raster(list.files(dir$source, full.names=T, recursive=T)[1])
mask <- crop(mask,r)
id <- mask
id[] <- 1:ncell(mask)
rid <- rotate(id)
rmask <- rotate(mask)

3. Creating Transformation Matrices

# 1. Extract overlap between raster and country polygon map (3 min) 
system.time( v <- extract(rid, map, small=TRUE, cellnumbers=TRUE, weights=TRUE))


# 2. Remove grid cells over water and rescale weights
v2 <- lapply(v, function(df) {
if (F){ # plot intermediate data
i <- 1 # Algeria
df <- v[[i]]
temp <- rmask
temp[df[,'cell']] <- 2
plot(temp)    
table(rmask[df[,'cell']])
(map@data$ADMIN.1)[i]
}
df <- df[rmask[df[,'cell']]==1, , drop=F] # keep cells over land
df[,'weight'] <- df[,'weight'] / sum(df[,'weight']) # rescale to 1
df # return
})


# 3. Remove countries too small to have any grid cells

# Check which countries are these
i <- which(sapply(v2,nrow)==0)
cbind(as.character(map@data$NAME[i]))
plot(map[i,], col='red', border='red', lwd=5)
plot(map, add=T, lwd=.5, col='grey')

# Remove them from weights list and map
v2 <- v2[-i]
map <- map[-i,]


# 4. Create list of transofrmation matrices + create figure

# Loop over aggregation criteria
fname <- paste(dir$figures,'/climate_aggregation_criteria.pdf',sep='')
#png(fname, width=700, height=length(misc$criteria)*250, pointsize=20)
cairo_pdf(fname, onefile=T, width=12, height=12, pointsize = 25)
par(mfrow=c(length(misc$criteria),1), mar=c(0,0,1,0), oma=c(0,6,0,0))
plist <- lapply(misc$criteria, function(crit) { # crit <- c('cropland','cor')[2]

print(crit)

# 1. Get weights into a raster layer

  # Import land cover shares
r2 <- readRDS(paste(dir$weights,'/',crit,'_share_Princeton.RDS',sep=''))

  # Crop and rotate to mask climate data size and orientation
r2 <- crop(r2, rid)
r <- mask
r[id[]] <- r2[rid[]]

  # View map
colors <- c('grey80',colorRampPalette(brewer.pal(9,'Reds'))(99))
plot(rotate(r), box=F, axes=F, legend=F, main=misc$criterianame[[crit]], col=colors)
plot(map, add=T, lwd=1, border='grey20')


# 2. Create empty sparse matrix

p <- Matrix(0,ncol=ncell(mask), nrow=length(unique(map@data$ISO3)), sparse=TRUE)
row.names(p) <- map@data$ISO3

# 3. Loop over countries (rows of sparse matrix)
for (i in row.names(p)) { 
  #print(i)
if (F) {
    # Check what happens when countries have no corpland weights (small islands, etc)
i <- 'ISL'
temp1 <- map[match(i,map@data$ISO3),]
temp2 <- map[which(map@data$GEO3==map@data$GEO3[map@data$ISO3==i]),]
plot(crop(rmask,temp2))
plot(temp1, add=T, border='red', lwd=5)
plot(crop(rotate(r),temp2))
plot(temp1, add=T, border='red', lwd=5)
}

  # a. Get the criteria weights
pos <- match(i, map@data$ISO3)
cells <- v2[[pos]][,'value']
weights <- r[cells]

  # b. Repare if ANY weights are missing OR all weights sum to zero
if (sum(is.na(weights))>0 | sum(weights,na.rm=T)==0) {
    # If ALL weights are missing OR weights sum to zero
if (sum(is.na(weights)|sum(weights==0,na.rm=T))==length(weights)) {
      weights <- v2[[pos]][,'weight']
    }
    # If SOME weights are missing
if (sum(is.na(weights))!=length(weights)) {
      keep <- which(!is.na(weights))
      cells <- cells[keep]
      weights <- weights[keep]
    }
  }

  # c. Make sure weights scale to 1
weights <- weights/sum(weights)

  # d. Set weights to sparse matrix
p[paste(i),cells] <- weights
}

# 4. Double check
if (F) {
  temp <- mask
  temp[] <- NA
  cells <- colSums(p)
  temp[] <- (cells>0)*1
  plot(rotate(temp)) # show cells that are represented in P matrix
}

# 5. Return
return(p)

})
names(plist) <- misc$criteria
dev.off()

# 5. Write list of transformation matrices to disk

fname <- paste(dir$weights,'/plist.RDS',sep='')
saveRDS(plist, fname)  
# The end

Explanation

  1. Setup and Environment: The script begins by setting up the workspace, loading required packages, defining directory paths, and defining miscellaneous parameters for data analysis.

  2. Data Preparation: The script loads raster data representing a land cover mask and a polygon dataset containing country boundaries. The script crops the land cover mask to match the extent of the climate data, and performs other data processing steps to prepare the data for creating transformation matrices.

  3. Creating Transformation Matrices: The script iterates through various criteria (e.g., cropland shares, pasture shares, etc.) to create transformation matrices for different weighting schemes. The code identifies overlapping cells between the raster data and country polygons, removes cells over water, and normalizes weights. For countries without any overlapping cells, alternative weighting schemes are applied. The script then creates sparse matrices that represent the transformation matrices, where each row corresponds to a country and each column corresponds to a raster cell.

  4. Output: The script saves the list of transformation matrices as a .RDS file, which can be used for subsequent analysis and data extraction.

In summary: This script provides a comprehensive framework for generating transformation matrices to extract country-level climate data based on various weighting schemes. These matrices can be used to calculate spatially-weighted averages of climate variables, representing country-level climate indicators.

Creating Transformation Matrices for Extracting Country-Level Climate Data

原文地址: http://www.cveoy.top/t/topic/il2m 著作权归作者所有。请勿转载和采集!

免费AI点我,无需注册和登录