1

The data I am using was not generated by me, however, I am stuck with using it anyways. My issue is, that in R, the CRS of the raster stack is missing or NA. If I open it with ArcGIS or Unidata (!the file is a NetCDF!) it gives me an odd description I am not aware of. ('char lambert_conformal_conic_2sp') Overall within the NetCDF are 2 different layers, storing (SWI_001 and SWI_010) pixel values I want to extract for every single day of the year.

I tried to check this webpage for further information about the 'char lambert_conformal_conic_2sp' projection, but it led to nowhere (https://mygeodata.cloud/cs2cs/).

Of course I also tried several different coding approaches, in order to find a way to assign a CRS to the file, similar to:

proj4string(r)=CRS("+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.81977 +y_0=2121415.69617 +ellps=WGS84 +units=m +no_defs ")

f <- projectRaster (r, '+proj=longlat + datum=WGS84')

However, I also tried to assign the same CRS of my file to the poylogon I am trying to extract the final data with, but since the CRS of the NetCDF file is NA to begin with, it doesnt work.

fieldnew <- sf::st_transform(field, crs = "***") #same projection as raster stack

The code I am working with at the moment is:

library(raster)
library(sf)
library(rgdal)
library(exactextractr)
library(ncdf4)
library(ggplot2)
library(reshape2)

files <- list.files('D:\\***\\***\\ASCAT_2018',pattern='*.nc',full.names=TRUE)

r <- stack(files, varname='SWI_001')

field <- read_sf(dsn = 'D:\\***\\***\\Studyarea_Polygon\\Polygon for R', layer = "testfield")

##plot extend of one file
plot(r[[1]], ext=10*extent(fieldnew)) ## check if polygon aligns with raster and if values make sense
plot(fieldnew, add=TRUE)  

## extract by weights of polygon Raster,SpatialPoylgon and relevant observations
output <- extract(r [[1:197]], fieldnew, fun=mean, weigths=TRUE, normalizeWeights=FALSE, layer=1, method='bilinear', df=TRUE)

The error messages are:

> f <- raster(files, varname='SWI_001')
Warning:
1: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:
long_name=CRS definition
spatial_ref=PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6377397.16,299.1527033239178]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",46],PARAMETER["standard_parallel_2",49],PARAMETER["latitude_of_origin",47.5],PARAMETER["central_meridian",13.333333],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
GeoTransform=-380501 500 0 220501 0 -500
2: In .getCRSfromGridMap4(atts) : cannot create a valid CRS
grid_mapping_name; false_easting; false_northing; scale_factor_at_projection_origin; scale_factor_at_central_meridian; standard_parallel; standard_parallel1; standard_parallel2; longitude_of_central_meridian; longitude_of_projection_origin; latitude_of_projection_origin; straight_vertical_longitude_from_pole; longitude_of_prime_meridian; semi_major_axis; semi_minor_axis; inverse_flattening; earth_radius; +proj; +x_0; +y_0; +k_0; +k_0; +lat_1; +lat_1; +lat_2; +lon_0; +lon_0; +lat_0; +lon_0; +pm; +a; +b; +rf; +a

Is there a way to assign a CRS to my files, or is it a lost cause?

EDIT: After futher research I found some entries about gdal_warp (gdalUtilis) and CDO (Changing a raster stack projection in R), is it possible to assign/transfrom the missing CRS based on these options?

3
  • Have you tried st_crs(name_of_your_raster)=4326? Commented Feb 4, 2021 at 14:45
  • @coconn41 thank you, I found the solution in QGIS. Checking the Metadata and name of the CRS, after I assigned it to my polygon and it worked for my process.
    – daves_ma
    Commented Feb 9, 2021 at 9:17
  • @coconn41, but I will try it with my data to have a new approach for issues like that.
    – daves_ma
    Commented Feb 9, 2021 at 9:18

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.