R Notebook

Ya se tiene descargado los datos de precipitación formato tif

De los global-pentad nombre: Chirps.v.2.0 2020.03.6

Arreglo de los datos

library(rgdal)

Loading required package: sp
rgdal: version: 1.5-10, (SVN revision 1006)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
Path to GDAL shared files: C:/Users/amor/Documents/R/win-library/4.0/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: C:/Users/amor/Documents/R/win-library/4.0/rgdal/proj
Linking to sp version:1.4-2
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

library(raster)
library(sf)

Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.2     v purrr   0.3.4
v tibble  3.0.1     v dplyr   1.0.0
v tidyr   1.1.0     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x tidyr::extract() masks raster::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x dplyr::select()  masks raster::select()

library(tmap)
library(gstat)
library(sp)

3. Preprocesamiento de datos CHIRPS

Lea el archivo CHIRPS sin comprimir:

precip <- raster("C:/Users/amor/Desktop/Rstudio/chirps-v2.0.2020.03.6.tif")

Como ya habrá notado, este archivo representa la lluvia acumulada en los últimos días de marzo de 2020 (es decir, la precipitación del 26 al 31 de marzo).

Verifique el contenido del objeto precip

precip

class      : RasterLayer 
dimensions : 2000, 7200, 14400000  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : C:/Users/amor/Desktop/Rstudio/chirps-v2.0.2020.03.6.tif 
names      : chirps.v2.0.2020.03.6 

Tenga en cuenta que este es un conjunto de datos con cobertura global. Su CRS es un sistema de coordenadas geográficas.

Carguemos un archivo shape que represente nuestra área de interés:

(aoi <- shapefile("C:/Users/amor/Desktop/Rstudio/MGN2017_50_META/50_META/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 29 
extent      : -74.89921, -71.07753, 1.604238, 4.899101  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO,   MPIO_CNMBR,                           MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,    Shape_Leng,       Shape_Area 
min values  :         50,      50001,     ACACÍAS,                                 1840,  117.34108343,      2017,       META, 0.43640580036, 0.00955216818621 
max values  :         50,      50711, VISTAHERMOSA, Ordenanza 63 de Noviembre 21 de 1965, 17247.6864753,      2017,       META, 8.63318071739,    1.40216640329 

Tenga en cuenta que los sistemas de referencia de coordenadas para ambos conjuntos de datos son los mismos.

Ahora, recortemos los datos de precipitación:

# Datos de precipitación de cultivos por extensión de área de interés
precip.crop <- raster::crop(precip, extent(aoi))

Ahora, enmascarar la capa de trama:

precip.mask <- mask(x = precip.crop, mask = aoi)

Verifique la salida:

precip.mask 

class      : RasterLayer 
dimensions : 66, 76, 5016  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -74.9, -71.1, 1.599999, 4.899999  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : chirps.v2.0.2020.03.6 
values     : 0.9892889, 34.25411  (min, max)

Luego, se plotea la capa Raster recortada:

plot(precip.mask, main= "CHIRPS precipitaciones en Meta desde 26.03. hasta 30.03 en 2020 [mm] ")
plot(aoi, add=TRUE)

Se puede hacer un mejor mapa con Leaflet usando la función addRasterImage:

library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precip.mask),
  na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(precip.mask, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(precip.mask),
    title = "CHIRPS precipitaciones en Meta desde 26.03. hasta 30.03 en 2020 [mm]")

Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 strings

La conversión de ráster en puntos se puede realizar utilizando la función rasterToPoints de la biblioteca ráster:

precip.points <- rasterToPoints(precip.mask, spatial = TRUE)

precip.points

class       : SpatialPointsDataFrame 
features    : 2765 
extent      : -74.875, -71.125, 1.624999, 4.824999  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       : chirps.v2.0.2020.03.6 
min values  :     0.989288926124573 
max values  :      34.2541122436523 

names(precip.points) <- "Lluvia"

precip.points

class       : SpatialPointsDataFrame 
features    : 2765 
extent      : -74.875, -71.125, 1.624999, 4.824999  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :            Lluvia 
min values  : 0.989288926124573 
max values  :  34.2541122436523 

str(precip.points)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 2765 obs. of  1 variable:
  .. ..$ Lluvia: num [1:2765] 8 7.8 6.26 7.17 8.09 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:2765, 1:2] -71.2 -71.1 -71.3 -71.2 -71.2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -74.87 1.62 -71.12 4.82
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"

Tracemos los puntos:

plot(precip.mask, main= "CHIRPS precipitaciones en Meta desde 26.03. hasta 30.03 en 2020 [mm]")
plot(aoi, add=TRUE)

points(precip.points$x, precip.points$y, col = "red", cex = .4)

Escribamos los objetos espaciales de puntos en un archivo de disco. Solo para ahorrar tiempo y trabajo en caso de que algo salga mal.

geojsonio::geojson_write(precip.points, file = "C:/Users/amor/Desktop/Rstudio/ppoints.geojson")

Success! File is at C:/Users/amor/Desktop/Rstudio/ppoints.geojson

<geojson-file>
  Path:       C:/Users/amor/Desktop/Rstudio/ppoints.geojson
  From class: SpatialPointsDataFrame

Leamos los puntos de precipitación. Consulte la documentación de geojsonio para comprender qué parámetros deben pasarse a la función geojson_read.

precip.points <- geojsonio::geojson_read("C:/Users/amor/Desktop/Rstudio/ppoints.geojson", what="sp")

Verifiquemos cuál es el objeto de puntos precip.:

precip.points

class       : SpatialPointsDataFrame 
features    : 2765 
extent      : -74.875, -71.125, 1.624999, 4.824999  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :            Lluvia 
min values  : 0.989288926124573 
max values  :  34.2541122436523 

¿Cuál es el punto de convertir los datos de precipitación ráster en datos de precipitación puntual?

Bueno, hay pocas estaciones meteorológicas de la Organización Meteorológica Mundial (OMM) en nuestro país:

Por lo tanto, CHIRPS puede ser una opción para obtener los datos de precipitación puntuales necesarios para obtener una superficie de precipitación continua.

Preparación adicional

En caso de que se pierda la conexión, podemos leer nuevamente el archivo de forma con el área de interés:

(aoi <- shapefile("C:/Users/amor/Desktop/Rstudio/MGN2017_50_META/50_META/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 29 
extent      : -74.89921, -71.07753, 1.604238, 4.899101  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO,   MPIO_CNMBR,                           MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,    Shape_Leng,       Shape_Area 
min values  :         50,      50001,     ACACÍAS,                                 1840,  117.34108343,      2017,       META, 0.43640580036, 0.00955216818621 
max values  :         50,      50711, VISTAHERMOSA, Ordenanza 63 de Noviembre 21 de 1965, 17247.6864753,      2017,       META, 8.63318071739,    1.40216640329 

Necesitamos convertirlo a una característica espacial:

meta_sf <-  sf::st_as_sf(aoi)

Ahora, disuelva los límites internos:

(border_sf <-
  meta_sf %>%
  summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.89921 ymin: 1.604238 xmax: -71.07753 ymax: 4.899101
geographic CRS: WGS 84
      area                       geometry
1 85524.79 POLYGON ((-73.44639 2.38595...

Convertir de característica espacial a dataframe espacial:

 (border <- as(border_sf, 'Spatial')) 

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -74.89921, -71.07753, 1.604238, 4.899101  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :           area 
value       : 85524.79262911 

otra conversion

(meta.sf <- st_as_sf(aoi) %>% mutate(MUNIC = MPIO_CNMBR, CODIGO = MPIO_CCDGO) %>% select(MUNIC, CODIGO))

Simple feature collection with 29 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.89921 ymin: 1.604238 xmax: -71.07753 ymax: 4.899101
geographic CRS: WGS 84
First 10 features:
               MUNIC CODIGO                       geometry
1      VILLAVICENCIO  50001 POLYGON ((-73.69288 4.28821...
2      ACACÃ\u008dAS  50006 POLYGON ((-73.80268 4.19734...
3   BARRANCA DE UPIA  50110 POLYGON ((-73.01607 4.63503...
4           CABUYARO  50124 POLYGON ((-73.08865 4.46048...
5  CASTILLA LA NUEVA  50150 POLYGON ((-73.7252 3.91777,...
6           CUBARRAL  50223 POLYGON ((-74.20143 4.0118,...
7            CUMARAL  50226 POLYGON ((-73.55407 4.31821...
8        EL CALVARIO  50245 POLYGON ((-73.74639 4.45008...
9        EL CASTILLO  50251 POLYGON ((-73.92004 3.74881...
10         EL DORADO  50270 POLYGON ((-73.83554 3.76362...

Ahora, intentaremos una tarea st_intersection. (https://cran.r-project.org/web/packages/sf/vignettes/sf3.html#geometrical_operations) se hace usando la biblioteca sf

# conversion
p.sf <- st_as_sf(precip.points)

## intersecta polígonos con puntos, manteniendo la información de ambos
(precip.sf = st_intersection(meta.sf, p.sf))

although coordinates are longitude/latitude, st_intersection assumes that they are planar
attribute variables are assumed to be spatially constant throughout all geometries

Simple feature collection with 2765 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -74.875 ymin: 1.624999 xmax: -71.125 ymax: 4.824999
geographic CRS: WGS 84
First 10 features:
              MUNIC CODIGO   Lluvia                 geometry
29   PUERTO GAITÃ\u0081N  50568 7.996710 POINT (-71.175 4.824999)
29.1 PUERTO GAITÃ\u0081N  50568 7.801209 POINT (-71.125 4.824999)
29.2 PUERTO GAITÃ\u0081N  50568 6.263782 POINT (-71.275 4.774999)
29.3 PUERTO GAITÃ\u0081N  50568 7.165887 POINT (-71.225 4.774999)
29.4 PUERTO GAITÃ\u0081N  50568 8.093927 POINT (-71.175 4.774999)
29.5 PUERTO GAITÃ\u0081N  50568 8.186864 POINT (-71.125 4.774999)
29.6 PUERTO GAITÃ\u0081N  50568 5.463953 POINT (-71.425 4.724999)
29.7 PUERTO GAITÃ\u0081N  50568 8.316946 POINT (-71.375 4.724999)
29.8 PUERTO GAITÃ\u0081N  50568 6.087829 POINT (-71.325 4.724999)
29.9 PUERTO GAITÃ\u0081N  50568 6.038989 POINT (-71.275 4.724999)

Dos tareas reproyecciones

p.sf.magna <- st_transform(precip.sf, crs=3116)

meta.sf.magna <- st_transform(meta.sf, crs=3116)

Una conversion

(precip2 <- as(p.sf.magna, 'Spatial'))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

class       : SpatialPointsDataFrame 
features    : 2765 
extent      : 911335.8, 1328392, 671451.4, 1026012  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 3
names       :        MUNIC, CODIGO,            Lluvia 
min values  :     ACACÍAS,  50001, 0.989288926124573 
max values  : VISTAHERMOSA,  50711,  34.2541122436523 

Nuevamente, podemos escribir el conjunto de datos intermedio (por si acaso):

shapefile(precip2, filename='C:/Users/amor/Desktop/Rstudio/precip2.shp', overwrite=TRUE)

Si es necesario, léelo. De lo contrario, sigue adelante

precip2$precipitaciones<- round(precip2$Lluvia, 1)

Lo que tenemos

precip2

class       : SpatialPointsDataFrame 
features    : 2765 
extent      : 911335.8, 1328392, 671451.4, 1026012  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :        MUNIC, CODIGO,            Lluvia, precipitaciones 
min values  :     ACACÍAS,  50001, 0.989288926124573,               1 
max values  : VISTAHERMOSA,  50711,  34.2541122436523,            34.3 

(meta2 <- as(meta.sf.magna, 'Spatial'))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

class       : SpatialPolygonsDataFrame 
features    : 29 
extent      : 908647.3, 1333690, 669155.7, 1034240  (xmin, xmax, ymin, ymax)

CRS object has comment, which is lost in output

crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :        MUNIC, CODIGO 
min values  :     ACACÍAS,  50001 
max values  : VISTAHERMOSA,  50711 

Es posible que desee escribir el nuevo objeto (por si acaso):

shapefile(meta2, filename='C:/Users/amor/Desktop/Rstudio/meta2.shp', overwrite=TRUE)

Asegúrese de que las dos extensiones coincidan:

# Replace point boundary extent with that of meta
precip2@bbox <- meta2@bbox

División de datos en conjuntos de datos de capacitación y validación

train_index <- sample(1:nrow(precip2), 0.75 * nrow(precip2))
test_index <- setdiff(1:nrow(precip2), train_index)
ptos_train <- precip2[train_index, ]
ptos_test  <- precip2[test_index,]
ptrain <- spTransform(ptos_train, crs(precip.mask))
ptest <- spTransform(ptos_test, crs(precip.mask))

#install.packages("htmltools")
library(htmltools)
library(leaflet.extras)

# Primero preparar leaflet plot ...
lplot <- leaflet(data = precip2) %>% # data = cuerpo original - para obtener el zoom correcto
  addProviderTiles("CartoDB.Positron") %>% 
  addRasterImage(precip.mask, colors = pal, opacity = 0.6) %>%
  addCircleMarkers(data = ptrain, # primer grupo
                   radius = 1,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(precipitaciones),
                   color = pal(ptos_train$precipitaciones), # usando una paleta ya creada
                   clusterOptions = markerClusterOptions(),
                   group = "Training") %>% 
  addCircleMarkers(data = ptest, # segundo grupo
                   radius = 10,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(precipitaciones),
                   color = pal(ptos_test$precipitaciones), # using already created palette
                   clusterOptions = markerClusterOptions(),
                   group = "Test") %>% 
  addLegend(position = "bottomright",
            values = ~precipitaciones,
            opacity = .7,
            pal = pal, # palette declared previously
            title = "CHIRPS precipitaciones en Meta \ndesde 26.03. hasta 30.03 en 2020 [mm]") %>% 
  leaflet::addLayersControl(overlayGroups = c("Training", "Test"),
                   options = layersControlOptions(collapsed = FALSE)) %>% 
  addResetMapButton()

Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defsDiscarded datum WGS_1984 in CRS definitionUsing PROJ not WKT2 stringsSome values were outside the color scale and will be treated as NASome values were outside the color scale and will be treated as NA

lplot #  ... display the map