search

r-spatial / stars

Spatiotemporal Arrays, Raster and Vector Data Cubes
https://r-spatial.github.io/stars/
Apache License 2.0
564 stars 94 forks source link

Setting ob_tran CRS #52

Closed adrfantini closed 5 years ago

adrfantini commented 6 years ago

I have a netcdf file (link) of which I know the proper CRS, but the CRS is not encoded in the file itself: +proj=ob_tran +o_proj=longlat +o_lat_p=29.0 +o_lon_p=0.0 +lon_0=15.0 +a=6367470 +b=6367470

stars, however, refuses to use it:

st_crs(s) = '+proj=ob_tran +o_proj=longlat +o_lat_p=29.0 +o_lon_p=0.0 +lon_0=15.0 +a=6367470 +b=6367470'
GDAL cannot import PROJ.4 string `+proj=ob_tran +o_proj=longlat +o_lat_p=29.0 +o_lon_p=0.0 +lon_0=15.0 +a=6367470 +b=6367470': returning missing CRS

This is very similar to what happened in r-spatial/sf#651 and has to do, as far as I understand it, with limited support in GDAL for the transformation ob_tran.

In the case of sf, a workaround using lwgeom was possible. Is there anything similar that can be done for stars?

edzer commented 6 years ago

Currently, we can read this file as a curvilinear grid, but when read it as a regular GDAL file with read_stars there is no information about the georeference, at all.

edzer commented 6 years ago

stars is not alone here:

> r = readGDAL("hmr_big.nc")
hmr_big.nc has GDAL driver netCDF 
and has 1361 rows and 1286 columns
Warning message:
In readGDAL("hmr_big.nc") : GeoTransform values not available

so in order to test this, please provide a netcdf file that at least gives a longlat grid, not just a matrix, when read through GDAL. Or give the key to this (offset, delta for x and y). As far as I understand, we shouldn't read these from the long and lat matrices, because they would contain the non-obtran transformed coordinates (ie the curvilinear grid).

adrfantini commented 6 years ago

Ok, I'll test this out and find a proper example in a few days.

adrfantini commented 6 years ago

Here you can find a time slice of an example netCDF file from one EURO-CORDEX model. The grid is a rotated one, defined as:

        int rotated_latitude_longitude ;
                rotated_latitude_longitude:grid_mapping_name = "rotated_latitude_longitude" ;
                rotated_latitude_longitude:grid_north_pole_latitude = 39.25 ;
                rotated_latitude_longitude:grid_north_pole_longitude = -162. ;
                rotated_latitude_longitude:north_pole_grid_longitude = 0. ;

lat and lon arrays are included. The CRS is not indicated explicitly. GDAL cannot understand the SRS of the file.

Here is another similar example.

Can this help?

edzer commented 6 years ago

What would be the proj4string for these datasets?

adrfantini commented 6 years ago

In theory it should be the same as previously discussed in r-spatial/sf#651, so something along the lines of "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329"

edzer commented 6 years ago

This uses the forward projection trick also used in the issue you refer to. 

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
(r = read_stars("example_ob_tran.nc"))
# stars object with 3 dimensions and 1 attribute
# attribute(s):
#  example_ob_tran.nc 
#  Min.   :247.0      
#  1st Qu.:267.4      
#  Median :275.9      
#  Mean   :274.5      
#  3rd Qu.:282.6      
#  Max.   :292.6      
# dimension(s):
#      from  to              offset delta  refsys point values    
# x       1 106              -28.43  0.44      NA    NA   NULL [x]
# y       1 103               21.89 -0.44      NA    NA   NULL [y]
# time    1   1 1979-01-16 12:00:00    NA POSIXct    NA   NULL

First the path using the lon and lat subdatasets to directly create a curvilinear grid:

lon = read_stars('NETCDF:"example_ob_tran2.nc":lon')
lat = read_stars('NETCDF:"example_ob_tran2.nc":lat')
rr = st_as_stars(r, curvilinear = list(x = lon[[1]], y = lat[[1]]))
rr
# stars object with 3 dimensions and 1 attribute
# attribute(s):
#  example_ob_tran.nc 
#  Min.   :247.0      
#  1st Qu.:267.4      
#  Median :275.9      
#  Mean   :274.5      
#  3rd Qu.:282.6      
#  Max.   :292.6      
# dimension(s):
#      from  to              offset delta                       refsys point
# x       1 106                  NA    NA +proj=longlat +datum=WGS8...    NA
# y       1 103                  NA    NA +proj=longlat +datum=WGS8...    NA
# time    1   1 1979-01-16 12:00:00    NA                      POSIXct    NA
#                                       values    
# x    [106x103] -44.1407 [°], ..., 64.404 [°] [x]
# y    [106x103] 22.1994 [°], ..., 72.4199 [°] [y]
# time                                    NULL    
# curvilinear grid
plot(rr, axes = TRUE, border = NA)

x1 Next compute the curvilinear grid parametrically from the shifted grid (r), using the forward trick:

obtran = "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +to_meter=0.01745329"
st_crs(r) = 4326
r.ll = st_transform(r, obtran)
st_crs(r.ll) = 4326
plot(r.ll, axes = TRUE, border = NA)

x2

r.ll
# stars object with 3 dimensions and 1 attribute
# attribute(s):
#  example_ob_tran.nc 
#  Min.   :247.0      
#  1st Qu.:267.4      
#  Median :275.9      
#  Mean   :274.5      
#  3rd Qu.:282.6      
#  Max.   :292.6      
# dimension(s):
#      from  to              offset delta                       refsys point
# x       1 106                  NA    NA +proj=longlat +datum=WGS8...    NA
# y       1 103                  NA    NA +proj=longlat +datum=WGS8...    NA
# time    1   1 1979-01-16 12:00:00    NA                      POSIXct    NA
#                               values    
# x    [106x103] -44.1407, ..., 64.404 [x]
# y      [106x103] 22.1994, ..., 72.42 [y]
# time                            NULL    
# curvilinear grid

st_bbox(rr)
# Units: [°]
#      xmin      ymin      xmax      ymax 
# -44.14069  22.19937  64.40398  72.41995 
st_bbox(r.ll)
#      xmin      ymin      xmax      ymax 
# -44.14070  22.19937  64.40399  72.41996
adrfantini commented 6 years ago

I guess that now after #65 the another approach would be: r2 = read_stars('example_ob_tran.nc', curvilinear=c('lon', 'lat')) Is one approach more 'correct' than the other?

edzer commented 6 years ago

the obtran PROJ string does not come from the file, lon and lat do. Using only stuff that comes from the file feels easier. Having said that, the obtran parameters again may be somewhere in the metadata. lon and lat arrays only work when they are present.

adrfantini commented 5 years ago

I am having issues with curvilinear for this file with rotated grid.

fn = 'historical.nc'
fn %>% read_stars(sub = 'Q100', curvilinear=c('lon', 'lat'))
#Q100, trying to read file: netCDF:NETCDF:"/home/clima-archive4-b/afantini/regcm_simulations/EURO-CORDEX/change/italy_change/mrro/faster_way/data_hadgem/yearmax/historical.nc":Q100:lon
#Error in CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  : 
  file not found

fn %>% read_ncdf(var = 'Q100', curvilinear=c('lon', 'lat'))
#Error: length(d) == length(dim(x[[1]])) is not TRUE

lon = fn %>% read_ncdf(var = 'lon', quiet = TRUE)
lat = fn %>% read_ncdf(var = 'lat', quiet = TRUE)
ll = setNames(c(lon, lat), c("x", "y"))
fn %>% read_ncdf(var = 'Q100', quiet = TRUE)  %>% st_as_stars(curvilinear = ll)
#Error: length(d) == length(dim(x[[1]])) is not TRUE

Can you reproduce this?

edzer commented 5 years ago 
fn = 'historical.nc'
lon = fn %>% read_ncdf(var = 'lon', quiet = TRUE)
lat = fn %>% read_ncdf(var = 'lat', quiet = TRUE)
lo = st_set_dimensions(lon, names = c("x", "y"))
la = st_set_dimensions(lat, names = c("x", "y"))
ll = setNames(c(lo, la), c("x", "y"))
x = fn %>% read_stars(sub = 'Q100', quiet = TRUE)  %>% st_as_stars(curvilinear = ll)

Any suggestions how we can make this easier?

And, any idea what went wrong? x

adrfantini commented 5 years ago

Shouldn't we be able to read this by just doing fn %>% read_stars(sub = 'Q100', curvilinear = c('lon', 'lat'))?

adrfantini commented 5 years ago

If you read lat and lon with read_stars instead of read_ncdf, it plots fine. Alternatively, I assume you could read everything with read_ncdf, but I am getting a Error: length(d) == length(dim(x[[1]])) is not TRUE if trying to plot or display x = fn %>% read_ncdf(var = 'Q100', quiet = TRUE) %>% st_as_stars(curvilinear = ll)

edzer commented 5 years ago

Thanks for sorting out! @mdsumner this may be related that in lon and lat both x and y dimensions have a negative delta (cellsize)?

mdsumner commented 5 years ago

Is that the GDAL default? I remember it would always be upside down if the "default" transform was applied, where offsets both zero and deltas both 1. It was the wrong orientation for a simple png image. I'll explore

edzer commented 5 years ago

First time I see a negative x delta! It is simply what the netcdf gdal driver communicates; I guess it comes from the ordering of the array in the netcdf file which is, afaict, not prescribed.

edzer commented 5 years ago 
fn = 'historical.nc'
fn %>% read_stars(sub = 'Q100', curvilinear = c('lon', 'lat')) %>% plot

now works.

adrfantini commented 5 years ago

Mh, not for me with GDAL 2.2.2. I am getting:

Q100, trying to read file: netCDF:NETCDF:"/dev/shm/afantini/historical.nc":Q100:lon
Error in CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  : 
  file not found
In addition: Warning message:
In CPL_read_gdal(as.character(x), as.character(options), as.character(driver),  :
  GDAL Error 1: Failed to parse NETCDF: prefix string into expected 2, 3 or 4 fields.

2.2.2 gdalinfo of the file:

Driver: netCDF/Network Common Data Format
Files: historical.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#absorption_emission_time_step_in_hours=18
  NC_GLOBAL#asselin_filter_nu=0.0625
  NC_GLOBAL#atmosphere_time_step_in_seconds=24
  NC_GLOBAL#boundary_high_nudge=12
  NC_GLOBAL#boundary_layer_scheme=2
  NC_GLOBAL#boundary_low_nudge=4
  NC_GLOBAL#boundary_medium_nudge=8
  NC_GLOBAL#boundary_nspgd=40
  NC_GLOBAL#boundary_nspgx=40
  NC_GLOBAL#c3s_disclaimer=This data has been produced in the context of the
 C3S_34b_Lot1 and Lot 2 projects (PRINCIPLES/CORDEX4CDS) as a data
 provider for the Climate Data Store within the Copernicus Climate Change
 Service (C3S - https://climate.copernicus.eu/). While abiding by the
 highest scientific and technical standards, ECMWF cannot warrant that
 any information provided by the C3S will be entirely free from errors or
 omissions or that such errors or omissions can or will be rectified
 entirely. This applies to data from projects that continue to be
 developed, but are made publicly available for the purpose of feedback
 and testing. Some data or metadata may have been created or structured
 in files or formats that are not error-free. ECMWF accepts no
 responsibility with regard to such problems incurred as a result of
 using this data (see http://climate.copernicus.eu/disclaimer-privacy for
 the full disclaimer)
  NC_GLOBAL#CDI=Climate Data Interface version 1.9.5 (http://mpimet.mpg.de/cdi)                                                                                                                                 
  NC_GLOBAL#CDO=Climate Data Operators version 1.9.5 (http://mpimet.mpg.de/cdo)                                                                                                                                 
  NC_GLOBAL#chemical_aerosol_scheme_activated=0                                                                                                                                                                 
  NC_GLOBAL#climatic_ozone_input_dataset=0                                                                                                                                                                      
  NC_GLOBAL#comment=RegCM CORDEX EURO-CORDEX_30 run                                                                                                                                                             
  NC_GLOBAL#contact=esp@ictp.it                                                                                                                                                                                 
  NC_GLOBAL#convection_time_step_in_seconds=240                                                                                                                                                                 
  NC_GLOBAL#convective_lwp_as_large_scale=1                                                                                                                                                                     
  NC_GLOBAL#Conventions=CF-1.4                                                                                                                                                                                  
  NC_GLOBAL#CORDEX_domain=EUR-11                                                                                                                                                                                
  NC_GLOBAL#creation_date=2018-08-10T12:05:06Z                                                                                                                                                                  
  NC_GLOBAL#cumulus_cloud_model=1                                                                                                                                                                               
  NC_GLOBAL#cumulus_convection_scheme_lnd=5                                                                                                                                                                     
  NC_GLOBAL#cumulus_convection_scheme_ocn=5                                                                                                                                                                     
  NC_GLOBAL#diffusion_hgt_factor=1
  NC_GLOBAL#diurnal_cycle_sst_scheme=0
  NC_GLOBAL#driving_experiment=MOHC-HadGEM2-ES, historical, r1i1p1
  NC_GLOBAL#driving_experiment_name=historical
  NC_GLOBAL#driving_model_ensemble_member=r1i1p1
  NC_GLOBAL#driving_model_id=MOHC-HadGEM2-ES
  NC_GLOBAL#dynamical_core=1
  NC_GLOBAL#experiment=EUR-11
  NC_GLOBAL#experiment_id=historical
  NC_GLOBAL#frequency=year
  NC_GLOBAL#grid_factor=0.7495652249943909
  NC_GLOBAL#grid_size_in_meters=12000
  NC_GLOBAL#history=2019-01-21 15:54:12 CET : created by create_Qx_regcm.R ; Mon Jan 21 15:19:58 2019: cdo yearmax data_hadgem/historical.nc data_hadgem/yearmax/historical.nc
Mon Jan 21 15:17:30 2019: cdo -L sellonlatbox,6.5,19,36,48 -selyear,1976/2005 /home/clima-archive4/CORDEX2/EUR-11/ICTP-RegCM4-6/MOHC-HadGEM2-ES/historical//mrro_EUR-11_MOHC-HadGEM2-ES_historical_r1i1p1_ICTP-RegCM4-6_v1_day_1970060112-2006010112.nc data_hadgem/historical.nc
Thu Jan 10 10:18:03 2019: cdo -O -z zip remap,/home/netapp-clima/users/jciarlo/scripts/postproc_cordex/sftlf_EUR-11_ICHEC-EC-EARTH_rcp45_r3i1p1_DMI-HIRHAM5_v1_fx.nc,../../../011grid_CORDEX_weights.nc mrro_EUR-11_MOHC-HadGEM2-ES_historical_r1i1p1_ICTP-RegCM4-6_v1_day_1970060112-1980010112.nc mrro_EUR-11_MOHC-HadGEM2-ES_historical_r1i1p1_ICTP-RegCM4-6_v1_day_1970060112-1980010112.nc_remap.nc
2017-10-31 16:25:30 : Created by RegCM RegCM Model program
  NC_GLOBAL#ICTP_version_note=Archived on model native grid
  NC_GLOBAL#institute_id=ICTP
  NC_GLOBAL#institution=International Centre for Theoretical Physics
  NC_GLOBAL#ipcc_scenario_code=HISTORICAL
  NC_GLOBAL#lake_model_activated=0
  NC_GLOBAL#landsurface_model=clm4.5
  NC_GLOBAL#large_scale_cloud_fraction_scheme=0
  NC_GLOBAL#lateral_boundary_condition_scheme=5
  NC_GLOBAL#latitude_of_projection_origin=48
  NC_GLOBAL#longitude_of_projection_origin=9.75
  NC_GLOBAL#model_icbc_data_source=HA_85
  NC_GLOBAL#model_id=ICTP-RegCM4-6
  NC_GLOBAL#model_is_restarted=No
  NC_GLOBAL#model_revision=tag-4.6.1
  NC_GLOBAL#model_simulation_end=1970-09-01 00:00:00 UTC
  NC_GLOBAL#model_simulation_initial_start=1970-06-01 00:00:00 UTC
  NC_GLOBAL#model_simulation_start=1970-06-01 00:00:00 UTC
  NC_GLOBAL#model_sst_data_source=HA_85
  NC_GLOBAL#moisture_scheme=1
  NC_GLOBAL#ncdf4_version=1.16
  NC_GLOBAL#note=The domain is larger than EUR-11
  NC_GLOBAL#ocean_flux_scheme=2
  NC_GLOBAL#pressure_gradient_scheme=0
  NC_GLOBAL#product=output
  NC_GLOBAL#projection=LAMCON
  NC_GLOBAL#project_id=CORDEX
  NC_GLOBAL#radiation_scheme_time_step_in_minuts=30
  NC_GLOBAL#rcm_version_id=v1
  NC_GLOBAL#references=http://gforge.ictp.it/gf/project/regcm
  NC_GLOBAL#rrtm_radiation_scheme_activated=0
  NC_GLOBAL#R_version=R version 3.4.4 (2018-03-15)
  NC_GLOBAL#seasonal_desert_albedo=0
  NC_GLOBAL#semi_lagrangian_advection_scheme=0
  NC_GLOBAL#simple_sea_ice_scheme=0
  NC_GLOBAL#source=RegCM Model output file
  NC_GLOBAL#standard_parallel={30,65}
  NC_GLOBAL#static_solar_constant_used=0
  NC_GLOBAL#subex_auto_conversion_rate_for_land=5e-05
  NC_GLOBAL#subex_auto_conversion_rate_for_ocean=0.00025
  NC_GLOBAL#subex_bottom_level_with_no_clouds=1
  NC_GLOBAL#subex_cloud_fraction_maximum=0.75
  NC_GLOBAL#subex_cloud_fraction_max_for_convection=1
  NC_GLOBAL#subex_cloud_liqwat_max_for_convection=0.0003
  NC_GLOBAL#subex_condensation_threshold=1
  NC_GLOBAL#subex_gultepe_factor_when_rain_for_land=0.4
  NC_GLOBAL#subex_gultepe_factor_when_rain_for_ocean=0.4
  NC_GLOBAL#subex_land_raindrop_accretion_rate=3
  NC_GLOBAL#subex_land_raindrop_evaporation_rate=0.01
  NC_GLOBAL#subex_limit_temperature=238
  NC_GLOBAL#subex_ocean_raindrop_accretion_rate=3
  NC_GLOBAL#subex_ocean_raindrop_evaporation_rate=0.001
  NC_GLOBAL#subex_rh_threshold_for_land=0.9
  NC_GLOBAL#subex_rh_threshold_for_ocean=0.9
  NC_GLOBAL#subex_rh_with_fcc_one=1.01
  NC_GLOBAL#surface_emissivity_factor_computed=0
  NC_GLOBAL#surface_interaction_time_step_in_seconds=600
  NC_GLOBAL#tiedtke_actual_scheme=4
  NC_GLOBAL#tiedtke_cape_adjustment_timescale=10800
  NC_GLOBAL#tiedtke_cloud_cover_evap_over_land=0.05
  NC_GLOBAL#tiedtke_cloud_cover_evap_over_ocean=0.05
  NC_GLOBAL#tiedtke_cloud_water_conv_over_land=0.008999999999999999
  NC_GLOBAL#tiedtke_cloud_water_conv_over_ocean=0.003
  NC_GLOBAL#tiedtke_coeff_evap_over_land=5.55e-05
  NC_GLOBAL#tiedtke_coeff_evap_over_ocean=5.55e-05
  NC_GLOBAL#tiedtke_critical_rh_over_land=0.5
  NC_GLOBAL#tiedtke_critical_rh_over_ocean=0.6
  NC_GLOBAL#tiedtke_cumulus_downdraft=Yes
  NC_GLOBAL#tiedtke_cumulus_friction=Yes
  NC_GLOBAL#tiedtke_detrainment_rate_deep=7.499999999999999e-05
  NC_GLOBAL#tiedtke_entrainment_rate_deep=0.00175
  NC_GLOBAL#tiedtke_entrainment_rate_downdraft=0.0003
  NC_GLOBAL#tiedtke_ke_dissipation=Yes
  NC_GLOBAL#tiedtke_midlevel=Yes
  NC_GLOBAL#tiedtke_penetrative=Yes
  NC_GLOBAL#tiedtke_prognostic_cloud=Yes
  NC_GLOBAL#tiedtke_shallow=Yes
  NC_GLOBAL#tiedtke_shallow_entrainment=2
  NC_GLOBAL#tiedtke_shallow_wstar_closure=No
  NC_GLOBAL#tiedtke_tracer_smooth_massflux=No
  NC_GLOBAL#tiedtke_tracer_transport=Yes
  NC_GLOBAL#title=ICTP Regional Climatic model V4
  NC_GLOBAL#tracking_id=db7822a6-9c84-11e8-b848-0894ef50a14e
  NC_GLOBAL#uwpbl_advection_scheme=0
  NC_GLOBAL#uwpbl_cloud_evap_entr_incr_efficiency=15
  NC_GLOBAL#uwpbl_czero=5.869
  NC_GLOBAL#uwpbl_eddy_LS_stable_PBL_scaling=1.5
  NC_GLOBAL#uwpbl_nuk=5
  NC_GLOBAL#zeng_ocean_roughness_formula=1
  NC_GLOBAL#zeng_ocean_roughness_method=1
  NC_GLOBAL#_NCProperties=version=1|netcdflibversion=4.6.1|hdf5libversion=1.10.1
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"historical.nc":time_bnds
  SUBDATASET_1_DESC=[30x2] time_bnds (64-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"historical.nc":lon
  SUBDATASET_2_DESC=[114x94] longitude (64-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"historical.nc":lat
  SUBDATASET_3_DESC=[114x94] latitude (64-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"historical.nc":mrro
  SUBDATASET_4_DESC=[30x114x94] runoff_flux (32-bit floating-point)
  SUBDATASET_5_NAME=NETCDF:"historical.nc":Q100
  SUBDATASET_5_DESC=[114x94] Q100 (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)
edzer commented 5 years ago

What happens if you add driver = NULL to the read_stars command?

adrfantini commented 5 years ago

Same error.

edzer commented 5 years ago

This is after installing the branch with

remotes::install_github("r-spatial/stars", "values")

right?

adrfantini commented 5 years ago

Dang! Sorry, I missed that. Yep, can confirm it's working if using that branch. I believe we can close this? The Error: length(d) == length(dim(x[[1]])) is not TRUE with read_ncdf persists but maybe that's material for another issue?

edzer commented 5 years ago

Why don't we get the full time series in the stars object?

edzer commented 5 years ago

Ah, sorry, that only concerns the mrro variable. x

adrfantini commented 5 years ago

Ah, see the damned 360-day calendar issue right there.

edzer commented 5 years ago

I noticed - time for its own temporal reference system. Did you mention there was an R package to convert back and forth between calendar time and 360-day/year model time?

mdsumner commented 5 years ago

@edzer where do you see a negative x? I see default geotransform on time, and typical +x/-y on the rest

# xoffset, xdelta, yskew, yres, xskew, ydelta (internal GDAL geotransform arrangement)
$time
[1] 0 1 0 0 0 1

$lon
[1] -9.51  0.11  0.00 -2.20  0.00 -0.11

$lat
[1] -9.51  0.11  0.00 -2.20  0.00 -0.11

$mrro
[1] -9.51  0.11  0.00 -2.20  0.00 -0.11

$Q100
[1] -9.51  0.11  0.00 -2.20  0.00 -0.11

in gdalinfo terms

Driver: netCDF/Network Common Data Format
Files: historical.nc
Size is 94, 114
Coordinate System is `'
Origin = (-9.509999923167690,-2.200000113090583)
Pixel Size = (0.109999998923271,-0.110000002700671)

But, I'd assume to ignore the geotransform in curvilinear case anyway - right? (GDAL may or may not override the default, depending on vicarious heuristics). Here using raster (divorced from how I tackle this in ocean model data - ignoring the extent/crs, but analogous) and using the quadmesh grid/base plot: 

f <- "historical.nc"
r <- raster::raster(f, varname = "mrro")
coords <- raster::brick(raster::raster(f, varname = "lon"), 
                        raster::raster(f, varname = "lat"))
library(quadmesh)
mesh_plot(r, coords = coords)  
maps::map(add = TRUE)

image

adrfantini commented 5 years ago

I noticed - time for its own temporal reference system. Did you mention there was an R package to convert back and forth between calendar time and 360-day/year model time?

Yup, PCICt!

edzer commented 5 years ago

@mdsumner I'm sorry, I mixed up the offset and delta colums!

I tried to recreate the figure you made but failed and got this: x are packages CRAN versions?

mdsumner commented 5 years ago

Oh! Yes, so easy to mix up - worldfile order, GDAL order, raster extent order ...

Yes both raster/quadmesh are CRAN versions.

Hmm, I'm seeing that error in the netcdf pathway - so I'll investigate that(!): 

x <- read_ncdf(fn, curvilinear = c("lon", "lat"))
x
#Error in dim.stars(x) : length(d) == length(dim(x[[1]])) is not TRUE
edzer commented 5 years ago

First shot at PCICt support:

library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.7.0, GDAL 2.3.2, PROJ 5.2.0
fn = 'historical.nc'
(xx = read_stars(fn, sub = 'mrro', curvilinear = c("lon", "lat")))
# lon, 
# lat, 
# mrro, 
# stars object with 3 dimensions and 1 attribute
# attribute(s):
#  historical.nc":mrro [kg/m^2/s]
#  Min.   :0                     
#  1st Qu.:0                     
#  Median :0                     
#  Mean   :0                     
#  3rd Qu.:0                     
#  Max.   :0                     
#  NA's   :173400                
# dimension(s):
#      from  to     offset    delta                       refsys point
# x       1  94         NA       NA +proj=longlat +datum=WGS8...    NA
# y       1 114         NA       NA +proj=longlat +datum=WGS8...    NA
# time    1  30 1976-07-01 360 days                        PCICt    NA
#                                    values    
# x    [94x114] 3.91253 [°],...,19.1685 [°] [x]
# y     [94x114] 35.3469 [°],...,48.495 [°] [y]
# time                                 NULL    
# curvilinear grid
plot(xx[,,,1:12])
# Warning messages:
# 1: In seq(from = x$offset + (x$from - 1 + where) * x$delta, by = x$delta,  :
#   Incompatible methods ("+.PCICt", "Ops.difftime") for "+"
# 2: In seq(from = x$offset + (x$from - 1 + where) * x$delta, by = x$delta,  :
#   Incompatible methods ("+.PCICt", "Ops.difftime") for "+"

Look at the labels: x

adrfantini commented 5 years ago

Wonderful! Should probably continue this in #29