Plotting NCAR MUSICA outputs using R

-

 

library(ncdf4)

library(colorout)

library(ggplot2)

library(cptcity)


fc <- list.files(path = "/glade/scratch/sibarra/SAAP/run",

                 pattern = "SAAP.cam.h0.2014",

                 full.names = T)

nc <- list.files(path = "/glade/scratch/sibarra/SAAP/run",

                 pattern = "SAAP.cam.h0.2014")


dates <- as.Date(nc, format = "SAAP.cam.h0.%Y-%m-%d-00000.nc")


nd <- paste0(strftime(dates, "%Y%m%d"), ".nc")


mi <- map_data('world', wrap=c(0,360))



# geo ####
lat <- ncvar_get(nc1, varid = "lat")
lon <- ncvar_get(nc1, varid = "lon")


#PM25 ####

x <- ncvar_get(nc1, varid = "PM25")

dtx <- data.table::data.table(
  lat = lat,
  lon = lon,
  x = x[, 32]*1e9
)

summary(dtx)

rc <- classInt::classIntervals(var = dtx$x,
                               n = 100, 
                               style = "quantile")$brks

ggplot(dtx, 
       aes(x = lon, 
           y = lat, 
           colour = x),
       alpha = 0.3) +
  geom_point() +
  scale_colour_gradientn(expression(PM[2.5]~mu*g*m^{-3}),
                         colours = cpt(rev = T), 
                         values = rc/max(rc)) +
  geom_polygon(data = mi, 
               aes(x = long,
                   y = lat,
                   group = group),
               fill = NA, 
               colour = "grey50") +
  coord_quickmap() +
  labs(x = NULL, y = NULL, 
       title = paste0("NCAR MUSICA South America ", nc[32])) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(text = element_text(size = 20),
        legend.key.height = unit(4, "lines"),
        plot.title = element_text(hjust = 0.5))



rc <- classInt::classIntervals(var = dtx[lon > 250 & lon < 350 &
                                           lat > -60 & lat < 15]$x,
                               n = 100, 
                               style = "sd")$brks

ggplot(dtx[lon > 250 & lon < 350 &
             lat > -60 & lat < 15], 
       aes(x = lon, 
           y = lat, 
           colour = x),
       alpha = 0.3) +
  geom_point(size = 3) +
  scale_colour_gradientn(expression(PM[2.5]~mu*g*m^{-3}),
                         colours = lucky(rev = T), 
                         values = rc/max(rc)) +
  geom_polygon(data = mi, 
               aes(x = long,
                   y = lat,
                   group = group),
               fill = NA, 
               colour = "grey50") +
  coord_quickmap(xlim = c(250, 350), 
                 ylim = c(-60, 15), 
                 expand = c(0,0)) +
  labs(x = NULL, y = NULL, 
       title = paste0("NCAR MUSICA South America ", nc[32])) +
  theme(text = element_text(size = 20),
        legend.key.height = unit(4, "lines"),
        plot.title = element_text(hjust = 0.5))





Comentarios

  1. MaryCarrollmiércoles, octubre 15, 2025

    The visualization method greatly switch free game clarifies complex atmospheric data. Providing sample scripts makes it much easier for readers to duplicate the analysis.

    ResponderEliminar

Publicar un comentario

Entradas populares de este blog

ERROR: Grib2 file or date problem, stopping in edition_num. Can't ungrib.exe SOLVED

-
Imagen
Well, it seems that this problem can be very common. Usually, you have a directory for downloading the meteorology from NCEP/UCAR. Probably, you will have the script for downloading the data in the same directory. If you link the grib files based on the name there won't be any problems ./link_grib.csh /path/GFS* However, it could be possible to link the whole directory. You might do this when you have a whole year of meteorological data. ./link_grib.csh /path/ If you have scripts for downloading the met data in the same directory, THEN you will have the message: ERROR: Grib2 file or date problem, stopping in edition_num To avoid this error, just move the script for downloading the met data elsewhere and link the directory again. Finally, do ungrib.exe UPDATE 2018-12/15 (I mean 15 of December of 2018) I had the same problem again, and with the same data! and this solution did not work till i did this: ./link_grib.csh /path/GFS Yes only the first letters, in ...
Leer más

Error If you REALLY want Grib2 output from WRF, modify the arch/Config_new.pl - WRF 3.7 in Ubuntu (Elementary)

-
Imagen
Well, i was installing WRF Chem in a server and also in my local machine, following this great tutorial, and when i did ./configure, i had the following error. I had all the libraries installed and environmental variables declared. If you REALLY want Grib2 output from WRF, modify the arch/Config_new.pl script. Right now you are not getting the Jasper lib, from the environment, compiled into WRF. I searched and searched and lots of dudes had the same problem. The solution was ridiculous, was change  one line in VRFV3/arch/Config_new.pl, the line 223, from  $I_really_want_to_output_grib2_from_WRF = "FALSE" ; to  $I_really_want_to_output_grib2_from_WRF = "TRUE" ; and ./configure in WRFV3 and, voilá The funny was, when i was searching for the answer, lots of people had the same problem. I hope it helps
Leer más

ARCGIS: Cálculo de pendiente entre puntos en el espacio xyz

-
Imagen
Se tienen puntos en el espacio con coordenadas xyz. Estas corresponden a coordenadas UTM huso 19, en metros. Pendiente = Distancia Vertical / Distancia Horizontal  Distancia Vertical = z2 - z1 Distancia Horizontal = d = [(y2-y1)^2 + (x2-x1)^2]^0.5 Pendiente = (z2 - z1) / [(y2-y1)^2 + (x2-x1)^2]^0.5 Pendiente porcentual = Pendiente *100 Entonces; x1 y1 z1 x2 y2 z2 pendiente% 344150.34 6302722 520.897 344911.2 6302394 530.5 1.159 344150.34 6302722 520.897 343774.7 6302378 521.858 0.188 345811.8 6302371 536 345893.9 6302547 533 -1.544 Sirve para calcular pendiente de caminos, ríos, etc :) En este caso lo hice para una modelación de flujo vehicular de Santiago Estraus , por lo que la diferencias de las alturas está según el sentido del flujo, desde la altura del nodo inicial al nodo final, por eso hacia el sector oriente se aprecia mayor pendie...
Leer más