mikejohnson51 / climateR

An R 📦 for getting point and gridded climate data by AOI
https://mikejohnson51.github.io/climateR/
MIT License
165 stars 40 forks source link

Question and Help: On How to Plot Monthly Data into Yearly Data In Tempreture change Map using getTerra climate #57

Closed Heed725 closed 2 years ago

Heed725 commented 2 years ago

Hello Mike; I have a Question and an Issue....Ive been trying to make a clipped levelplot of Tempreture(tmax) using getTerra function but ive been having 2 key issues

The first one is on turning Monthly raster to Yearly raster ......because it gives me error on this particular part tmax.1989.2018=list.files("B:/TREND MAP DATA/changemaptmax",pattern="*.tif",recursive=T,full.names=T) tmax.1989.2018.stack=stack(tmax.1989.2018) names(tmax.1989.2018.stack)=c("1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018") myPal <- RColorBrewer::brewer.pal('YlOrRd', n=4) myTheme <- rasterTheme(region = myPal) levelplot(tmax.1989.2018.stack, par.settings = myTheme)

levelplot(tmax.1989.2018.stack,layout=c(6,5),names.attr=c("1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018"), par.settings=myTheme)

It says error Unequal Number

and The second one is to change colors from green to heat colors and clipping into yearly data like the output Below

My full code

library(shapefiles) library(sf) library(rasterVis) library(terra) library(sp) library(chirps) library(tidyverse) library(climateR) library(AOI) library(remotes) library(raster) library(colorRamps) library(RColorBrewer) library(lattice)

Inserting Dodoma urban shape file

Dodoma_urban.shp<-shapefile("C:/Users/Charles/Desktop/Dodoma season 02/R Madini/data/dodoma urban/ddm/dodoma_urban.shp") plot(Dodoma_urban.shp)

Creating simple feature

Dodoma_urban.shp.sf<- as(Dodoma_urban.shp, "sf") plot(Dodoma_urban.shp.sf, max.plot=10)

DODOMA_CSV<-read.csv("C:/Users/Charles/Desktop/Dodoma season 02/R Madini/data/dodoma urban/ddm/dodoma.csv")

set coordinates

coordinates(DODOMA_CSV)=~Long+Lat crs(DODOMA_CSV)="+proj=longlat +datum=WGS84 +no_defs" DODOMA_CSV.sf<- as(DODOMA_CSV, "sf")

Downloading TerraClimate data using getTerraClim function

Dodoma_tdata<-getTerraClim(Dodoma_urban.shp.sf, param="tmax", startDate="1998-01-01", endDate="2019-12-31")

stacking

Dodoma_tdata<-stack(Dodoma_tdata)

plot(Dodoma_tdata)

dates = seq.Date(as.Date('2018-01-01'),as.Date('2018-12-31'),by = 'month') nm = paste0('tmax',dates) names(Dodoma_tdata1) = nm writeRaster(Dodoma_tdata1,paste0('A:/Dodoma project 4yr/LST-MOD11A2/Tmax2018',filename = names(Dodoma_tdata1)),format = 'GTiff',bylayer = TRUE, overwrite = TRUE)

Dodoma_urban.tmax= mask(crop(Dodoma_tdata,Dodoma_urban.shp),Dodoma_urban.shp) plot(Dodoma_urban.tmax)

tmax

tmax.1989.2018=list.files("B:/TREND MAP DATA/changemaptmax",pattern="*.tif",recursive=T,full.names=T) tmax.1989.2018.stack=stack(tmax.1989.2018) names(tmax.1989.2018.stack)=c("1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018") myPal <- RColorBrewer::brewer.pal('YlOrRd', n=4) myTheme <- rasterTheme(region = myPal) levelplot(tmax.1989.2018.stack, par.settings = myTheme)

levelplot(tmax.1989.2018.stack,layout=c(6,5),names.attr=c("1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018"), par.settings=myTheme)

The output I want to make photo_2022-06-09_12-07-24

Heed725 commented 2 years ago

Help

mikejohnson51 commented 2 years ago

Hi @Heed725,

I believe this is what you're after. I don't have access to your files so I choose an AOI of Colorado, USA, and a shorted timeframe for a faster result.

You should be able to modify this too replicate that image.

Thanks!

library(terra)
library(rasterVis)
library(AOI)
library(climateR)
library(RColorBrewer)

# Get 6 years of tmax data for Colorado
test_data <- getTerraClim(
  AOI = aoi_get(state = "CO"),
  param = "tmax",
  startDate = "2014-01-01",
  endDate = "2019-12-31"
)
#> Spherical geometry (s2) switched off
#> Spherical geometry (s2) switched on

# We used 6 years, and want a annual max, so we need to apply "max" to every set of 12 layers
index = rep(1:6, each = 12)

annualMax = tapp(rast(test_data$terraclim_tmax),
                 index = index,
                 fun = "max")

# Plot
levelplot(annualMax,
          par.settings = rasterTheme(region = brewer.pal('YlOrRd', n = 9)),
          names.attr = as.character(c(2014:2019)))

Created on 2022-06-09 by the reprex package (v2.0.1)

Heed725 commented 2 years ago

First of Thank you so much for that Code ,One minor problem .....how about croping according to area of Interest how can you able to do that ?

Heed725 commented 2 years ago

Like let say you have a colorado shape file and cropping it to the output got.......like Im asking how can you be able to reach to the other end

mikejohnson51 commented 2 years ago

I am not entirely sure what you mean. climateR functions will crop to the AOI provided. If you want to further crop the returned results I would use terra::crop(input, AOI)

Heed725 commented 2 years ago

No I mean like Colorado boundary or New york Boundary ......and then cropping it according to its boundary not like box

Heed725 commented 2 years ago

*and then clipping

mikejohnson51 commented 2 years ago

ooooh the operation you are looking for is called "masking":

library(terra)
library(rasterVis)
library(AOI)
library(climateR)
library(RColorBrewer)

# Get 6 years of tmax data for New York
AOI = aoi_get(state = "NY")

test_data <- getTerraClim(
  AOI = AOI,
  param = "tmax",
  startDate = "2014-01-01",
  endDate = "2019-12-31"
)

# We used 6 years, and want a annual max, so we need to apply "max" to every set of 12 layers
index = rep(1:6, each = 12)

annualMax = tapp(rast(test_data$terraclim_tmax),
                 index = index,
                 fun = "max") |> 
  mask(vect(AOI))

# Plot
levelplot(annualMax,
          par.settings = rasterTheme(region = brewer.pal('YlOrRd', n = 9)),
          names.attr = as.character(c(2014:2019)))

Created on 2022-06-09 by the reprex package (v2.0.1)

Heed725 commented 2 years ago

Brilliant Question Two Last questions The first one is there a way to Mask it with a Shape file?

secondly; This AOI package out of USA can It do like A certain Region which is Admin Level 2 in other country let say Kenya or its Limited in Admin Level 0 in other country?

Rapsodia86 commented 2 years ago

Use library(sf) to read your shapefile:

library(terra) library(rasterVis) library(sf) library(climateR) library(RColorBrewer)

Get 6 years of tmax data for New York

KGZ <- st_read("X:/YOURPATH/gadm40_KGZ_0.shp")

test_data <- getTerraClim( KGZ, param = "tmax", startDate = "2014-01-01", endDate = "2019-12-31" )

We used 6 years, and want a annual max, so we need to apply "max" to every set of 12 layers

index = rep(1:6, each = 12)

annualMax = tapp(rast(test_data$terraclim_tmax), index = index, fun = "max") |> mask(vect("X:/YOURPATH/gadm40_KGZ_0.shp"))

Plot

levelplot(annualMax, par.settings = rasterTheme(region = brewer.pal('YlOrRd', n = 9)), names.attr = as.character(c(2014:2019))) image

Heed725 commented 2 years ago

One Minor thing In this Part on Masking It Gives me This error

annualMax = tapp(rast(test_data$terraclim_tmax), index = index, fun = "max") |> mask(vect(AOI))

Error in h(simpleError(msg, call)) : error in evaluating the argument 'mask' in selecting a method for function 'mask': error in evaluating the argument 'x' in selecting a method for function 'vect': object 'AOI' not found

Heed725 commented 2 years ago

Rplot

The Output is There.....The only thing i want is to mask It According to shape file given

Rapsodia86 commented 2 years ago

I just updated my post. Here, should be instead AOI: mask(vect(path)) or mask(vect(variable)l) Either a path to your shapefile or a variable name which is used for the read shapefile

Heed725 commented 2 years ago

Thank You so Much Mike and Raspodia Ive managed to reach Where I wanted teemp d

mikejohnson51 commented 2 years ago

Thanks for jumping in @Rapsodia86 ! Im glad it had a successful resolution :)