Home > Uncategorized > Coastal

Coastal

UPDATE:

cce asked for a chart of coastal locations broken out by their SST component and their Land componement.

That’s a little tricky, but once you figure out the mask logic its just addition and such.There are TWO masks: One mask contains ALL zeros, except for coast cells which are fraction of LAND.The other mask contains all zeros except the same coast cell are fraction of water. NOTE. the values here are the contribution of warming from coastal cells. Its SMALL. remember everything is area weighted. But you could look at the trends.

Red is the “land” portion, blue is the sea portion. orange is the sum. index is months from 1900.

UPdate: if you want to follow along after this I suggest you get hooked up to raster on R Forge. I’ll be testing out new capability. The code to date uses the CRAN release 1.4.10, but I need to move up to 1.5.2 to test

the new capability. However, baselining the old code probably needs to happen first. A few more bits to test and then we will be shifting to Forge releases of raster ( Cran requires formal testing, so be patient)

2-Sep-2010,  version 1.5-2

improved speed of gridDistance (JvE)
rasterToPoints can now take a matrix of values and return a brick with layer for each column (suggested by Steven Mosher)
writeRaser (internal saveAs function) takes more care as not to overwrite its own source file, even if overwrite=TRUE

A weird little benefit of doing things in raster. I can get a look at coastal

if we do

Coast <- Land+SST

the result will be ONLY those areas where they overlap

coast image

Thats because SST has NA where there is 100% Land, and land has NA where there is 100% Ocean

Add those togther and you get NA everywhere BUT the coast

Land<-getMask(file=LandWaterMask_File)

Inv<-getGhcnInventory()

# this line goes away

Inv<-as.GriddedInventory(Inv,Land)

Anom<-loadAnomalies(getwd())

Data<-intersect.InvAnomalies(Inv,Anom)

Anom<-Data$Anomalies

Inv<-Data$GridInventory

#####################################################

# this code becomes a one liner with the new ‘pointstoRaster’ basically

# we take the Anomalies which are station based, and aggregate them into a

# cell structure taking the average of all the stations (points) that map to the raster cell.

# see R Forge if you dont understand. stations are points(x,y) you aggregate those points

# into a raster (cells of lat/lon) using a “mean” function. so the work of as.Grid goes away and most

# of the lines that follow. We basically read a V2anomaly structure into a brick. slick.

CellMatrix<-as.GridCells(Anom,Inv)

m<-matrix(nrow=ncell(Land),ncol=ncol(CellMatrix))

b<-brick(Land)

clearValues(b)

m[as.numeric(row.names(CellMatrix)),]<-as.matrix(CellMatrix)

b<-setValues(b,m)

##############################

Weights <- area(b,na.rm=T,weight=T)

Temps<-Weights*b*Land

SST <- loadSST(getwd())

OceanPercentMask <- 1 – Land

Weights <- area(SST,na.rm=T,weight=T)

Ocean <- SST*Weights*OceanPercentMask

Coastal <-Ocean+Temps

Categories: Uncategorized
  1. cce
    September 6, 2010 at 8:04 PM

    Can you generate two time series from the overlapping gridcells? One using the land anomalies, and the other using the ocean?

  2. Steven Mosher
    September 6, 2010 at 8:30 PM

    Ya sure, I’ll do that tonight. Its fairly straight forward.

    • cce
      September 7, 2010 at 5:00 AM

      Thanks Mosher. Could you upload those time series somewhere? This is sort of what McIntyre is doing in his recent Hawaii, except in this case it is for all ocean/land areas. I’m looking for step changes in SST. I think it would be better if the land and ocean components weren’t weighted by area within their gridcells. That is, we assume that SST readings within the cell are a good measure/proxy of the land within the same cell and vice versa. Then we compare the two series, take their differenc, etc.

  3. Steven Mosher
    September 7, 2010 at 9:26 AM

    Ya sure, But we just have to clear on what the represent. Lets just go through the math

    You have a SST series. every month the grids reporting vary.
    We take every month and do the following.

    we take the cells reporting and we calculate the area of the cells reporting and divide that by the total of all cells. that is the cell weights. it sums to 1. Next we multiple the temps by the weight, then we multiply by the fraction of water in each cell. With the coast map, ALL water cells that are
    100% get 0%. so we are just multiplying the coast cells by the weight.
    so a coast cell will be the temp*fractional area* actual water area. I could also do the weight just as the weight of all coast cells. in that case the weights of the coasts will sum to 1.

    We do the same thing for the Land.

    There are of course cases where the land may not report its portion of the cell. or vice versa.

    I think for what you want to do, I should probably limit the analysis to cells where we have BOTH reports on a given month. Anyways the files should be in the drop box to the right.

  4. cce
    September 7, 2010 at 11:11 AM

    Yes, we’d need to keep it just to cells with both land and ocean that month. Easier said than done, but the way I think of it, you treat all the cells with water and land as “the world”. Find the anomalies based only on the ocean data. Find the average anomaly (i.e. treat them all as equal weight, regardless of what portion is ocean). Then do the same thing except with land.

  5. Steven Mosher
    September 7, 2010 at 11:50 AM

    Ok let me think. When I add landtemps to water I get a result where ONLY the cells that both report are present. I can take that result turn it into a it into a mask and reapply it

    Coast<-Land+Ocean # add cells if the both report, otherwsie NA
    coastmask<-is.logical(Coast) # true if the both report
    Land<-Land*coastmask # get you land temp ONLY those coast cells where they both report
    Ocean <- Ocean*coastmask
    Land <- area(Land,weight=T)*Land
    Ocean<-area(Ocean,weight-T)*Ocean

    And then Land and ocean will hold ONLY those coast cells where they both have valuesscaled by area of the cell

    and you have want you want with these two lines

    x<-cellStats(Land,sum)
    y<-cellStats(Ocean,sum)

    wa la. the beauty of a good API. brick math.

  1. No trackbacks yet.

Leave a comment