polymask.R 4.52 KB
Newer Older
Burkhardt Rockel's avatar
Burkhardt Rockel committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
polymask<-function(Cfile, Pfile, Mfile=NULL, multiply=1, skiplines=0, undef=-1.E20, ...) {

#
#  return value of the funcion is the 2D-mask
#

# read longitude and latitude data 
nc <- nc_open(Cfile)
if ( is.null(nc$var$rotated_pole$id) ) {
  # normal geographical grid
  mlon <- nc$dim$lon$vals
  mlat <- nc$dim$lat$vals
} else {
  # rotated grid
    # read in the coordinates of the rotated pole
    pollon  <- as.numeric(ncatt_get(nc, "rotated_pole", "grid_north_pole_longitude")$value)
    pollat  <- as.numeric(ncatt_get(nc, "rotated_pole", "grid_north_pole_latitude")$value)
    polgam  <- ncatt_get(nc, "rotated_pole", "north_pole_grid_longitude")
    if (polgam$hasatt){
        polgam <- polgam$value
    } else {
        polgam <- 0
    }
    mlon    <- nc$dim$rlon$vals
    mlat    <- nc$dim$rlat$vals

    lon <- ncvar_get(nc, "lon")
    lat <- ncvar_get(nc, "lat")

}
nc_close(nc)

# read polygon data
tmp <- read.table(Pfile, skip=skiplines)
Burkhardt Rockel's avatar
Burkhardt Rockel committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

# find end points of closed polygons
fnd999 <- which(tmp[,1] == -999)
if (length(fnd999)==0) { 
   print("Single Polygon file") 
   poly.ends<-which(duplicated(tmp))
   poly.starts<-poly.ends - poly.ends[1] + 1
   poly.number <- 1
} else {
   poly.ends <- c(fnd999 - 1,nrow(tmp))
   poly.starts <- c(1,fnd999 + 1)
   poly.number <- length(poly.ends)
}

# loop over number of closed polygons
for (npoly in 1:poly.number) {

poly.start<-poly.starts[npoly]
poly.end<-poly.ends[npoly]

plon<-tmp[,1][poly.start:poly.end]
plat<-tmp[,2][poly.start:poly.end]

Burkhardt Rockel's avatar
Burkhardt Rockel committed
58 59 60 61 62 63 64 65 66 67 68
plon[length(plon)+1] <- plon[1]
plat[length(plat)+1] <- plat[1]

# transform polygon data onto rotated grid, if necessary
if ( ! is.null(nc$var$rotated_pole$id) ) {
  poly <- geo2rot(pollon, pollat, plon, plat, polgam)
  plon <- poly$x
  plat <- poly$y
}

# determine the polygon area
Burkhardt Rockel's avatar
Burkhardt Rockel committed
69 70 71 72 73 74 75 76 77 78 79 80
  outgrid.tmp <- frac_in_polygon(mlon, mlat, plon, plat, multiply = multiply)
  if (npoly == 1 ) {
    outgrid <- outgrid.tmp
  }
  else
  {
    outgrid <- outgrid + outgrid.tmp
  }
  
}

outgrid[outgrid > 1] <-1
Burkhardt Rockel's avatar
Burkhardt Rockel committed
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
outgrid[outgrid==0] <- undef

if ( ! is.null(Mfile) ) {

if ( is.null(nc$var$rotated_pole$id) ) {
    mlondim <- ncdim_def("lon", "degrees_east", mlon, longname="longitude")
    mlatdim <- ncdim_def("lat", "degrees_north", mlat, longname="latitude")
    var <- ncvar_def("MASK", "0-1", list(mlondim, mlatdim), undef,longname="mask",prec="float")
    ncid <- nc_create(Mfile, var)
    ncatt_put(ncid,"lon","standard_name","longitude",prec="text")
    ncatt_put(ncid,"lon","standard_name","latitude",prec="text")
    ncatt_put(ncid,var,"_FillValue",undef,prec="float")
    ncatt_put(ncid,0,"Conventions","CF-1.0",prec="text")
    ncvar_put(ncid, var, outgrid)
    nc_close(ncid)
#    system(paste('ncwa -h -O -a rotated_pole ',Mfile,' ',Mfile,sep="")) # to make rotated_pole a scalar

} else {

    mlondim <- ncdim_def("rlon", "degrees", mlon,longname="longitude in rotated pole grid")
    mlatdim <- ncdim_def("rlat", "degrees", mlat,longname="latitude in rotated pole grid")
    poledim<- ncdim_def("rotated_pole","", 1:1, create_dimvar=FALSE )
    londef <- ncvar_def("lon", "degrees_east",list(mlondim,mlatdim),longname="longitude",prec="float")
    latdef <- ncvar_def("lat", "degrees_north",list(mlondim,mlatdim),longname="latitude",prec="float")
    gmap <- ncvar_def("rotated_pole", "", poledim, prec="char")
    var <- ncvar_def("MASK", "0-1", list(mlondim, mlatdim), undef,longname="mask",prec="float")
    ncid <- nc_create(Mfile, list(londef,latdef,gmap,var))
    ncatt_put(ncid,"rlon","standard_name","grid_longitude",prec="text")
    ncatt_put(ncid,"rlat","standard_name","grid_latitude",prec="text")
    ncatt_put(ncid,londef,"standard_name","longitude",prec="text")
    ncatt_put(ncid,latdef,"standard_name","latitude",prec="text")
    ncatt_put(ncid,gmap,"grid_mapping_name","rotated_latitude_longitude",prec="text")
    ncatt_put(ncid,gmap,"grid_north_pole_longitude",pollon,prec="float")
    ncatt_put(ncid,gmap,"grid_north_pole_latitude",pollat,prec="float")
    if (polgam != 0.) {
      ncatt_put(ncid,gmap,"north_pole_grid_longitude",polgam,prec="float")
    }
    ncatt_put(ncid,var,"coordinates","lon lat",prec="text")
    ncatt_put(ncid,var,"grid_mapping","rotated_pole",prec="text")
    ncatt_put(ncid,var,"_FillValue",undef,prec="float")
    ncatt_put(ncid,0,"Conventions","CF-1.0",prec="text")
    ncvar_put(ncid, londef, lon)
    ncvar_put(ncid, latdef, lat)
    ncvar_put(ncid, var, outgrid)
    nc_close(ncid)
#    system(paste('ncwa -h -O --no_cll_mth -a rotated_pole ',Mfile,' ',Mfile,sep="")) # to make rotated_pole a scalar
}

}

Burkhardt Rockel's avatar
Burkhardt Rockel committed
131
  return("Everything ok")
Burkhardt Rockel's avatar
Burkhardt Rockel committed
132 133 134

}  # end of function polymask