crs required vs optional?

[shospital]

I was checking the conversion functions and the values from https://polarwatch.noaa.gov/erddap/griddap/nsidcCDRice_nh_grid.htmlTable?longitude%5B(5837500.0):1:(-5337500.0)%5D%5B(-3837500.0):1:(3737500.0)%5D,latitude%5B(5837500.0):1:(-5337500.0)%5D%5B(-3837500.0):1:(3737500.0)%5D

since crs was optional and didn't throw errors, I expected that crs info was read correctly from myInfo but the values were incorrect.

Only when I added crs = "EPSG:3413", the values look close to those of the nsidc grid lat/lon, x and y values.

Could we either throw errors for no crs, or require crs to be specified?

[rmendels]

@shospital Yep, absolutely correct. As I said in my email the code is missing a lot of checks for bad or missing input. This iteration just wanted to get code that basically worked. Also not totally positive the code works with the like nine different ways that the projection is encoded - some give CRS number, some given PROJ string, some give WKT. That is some more testing I need to do. Glad it worked if you gave the CRS. But please do keep posting problems you find so I can track them and test for them. Just that after working on code for a long while I need a break from it or else I get in a mental rut and see what hope is there rather than what is there. Also I am trying to not force the user to specify the CS because for many users they would have no idea how to find that info.

[shospital]

@rmendels when users need to convert from or to lat lon, they would expect they are working with projections so it may be helpful to provide warning they need to specify crs code (or other forms). otherwise they may be working with wrong values without knowing. Perhaps more examples could be added to vignette?

I tried EPSG code and proj4text and both worked.

latlon_to_xy(myInfo, latitude, longitude, crs ="+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs")

latlon_to_xy(myInfo, latitude, longitude, crs ="EPSG:3413")

[rmendels]

@shospital yes I need to ad more checking to the code. In the meantime:

1. When I run it with no defined crs it does throw an error, not certain what you are seeing, Please give R code that reproduces what you see.
2. The particular dataset you are using has the projection under a name that I didn't have in the list (there are now 9 different names). That has been added to the code. Please update conversion.R to get the changes. It just worked for me for datasetID 'nsidcCDRice_nh_grid' that you reference above.

Let me know what happens. Most importantly we need to identify all the names that are used to store projection data. Right now I have:

  proj_strings <- c('proj4string', 'proj_crs_code', 'proj4text', 'projection',
                    'grid_mapping_epsg_code','grid_mapping_proj4',
                    'grid_mapping_proj4_params', 'grid_mapping_proj4text',
                    'WKT')

[shospital]

@rmendels

1. Here's what I ran

library(lubridate)
library(ncdf4)
library(rerddap)
library(rerddapUtils)
myURL <- 'https://polarwatch.noaa.gov/erddap/'
myInfo <- rerddap::info('nsidcG02202v4sh1day', url = myURL)
latitude <- c(31.1026717524309,31.1994138236332 )
longitude <- c(168.320422464133, 168.148751056332)

latlon_to_xy(myInfo, latitude, longitude) 
#           X         Y
# [1,] 4405883 -21313460
# [2,] 4478507 -21342022

latlon_to_xy(myInfo, latitude, longitude, crs="EPSG:3413")
#             X       Y
# [1,] -3837418 5837376
# [2,] -3812419 5837376

2. srid (spatial reference identifier), authority:code (EPSG:3413), crs are what I've seen so far.
   From this article: https://r.geocompx.org/reproj-geo-data, many are not encouraged any more, and the standards are WKT authority code (epsg:3413).

For polarwatch, we try to put as many attribtues as possible for users with diff software preference but I am learning that the following should be the required fields: WKT, 'proj4string', 'proj_crs_code', srid in erddap metadata. I reached out to Michael and others in projection world.

[rmendels]

@shospital the reason for the discrepancy is because the information in the file says it is not EPSG:3413. (it is "EPSG:3412") If a crs is given I don't check if it is what is in the file, will have to think about that. But the code returned the correct values for the given crs, as well as for the crs in the file. You can see the file crs as so:

myURL <- 'https://polarwatch.noaa.gov/erddap/'
myInfo <- rerddap::info('nsidcG02202v4sh1day', url = myURL)
proj_strings <- c('proj4string', 'proj_crs_code', 'proj4text', 'projection',
                    'grid_mapping_epsg_code','grid_mapping_proj4',
                    'grid_mapping_proj4_params', 'grid_mapping_proj4text',
                    'WKT')
crs_test <- intersect( proj_strings,  myInfo$alldata$NC_GLOBAL$attribute_name)
proj_crs_code_index <- which(myInfo$alldata$NC_GLOBAL$attribute_name == crs_test[1] )
myInfo$alldata$NC_GLOBAL$value[proj_crs_code_index]
"EPSG:3412"

If that is a typo in the file that is the kind of thing it is hard for the program to know - and you should let the people generating the files know. I am aware of the recommendation, I am also aware of what is actually in the files and until whomever goes to a standard way of encoding the crs info I need to go with what is there to be useful. I should add that the file also has the attribute "grid_mapping_proj4text" which produces "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs" and using that string the result is the same as for "EPSG:3412".

[shospital]

@rmendels 😰 you are right. I put the wrong id. i meant to add the north hemisphere. case resolved. 💥

[rmendels]

@shospital there is new version of conversion.R. Please update to that and test. The logic has been refactored and there are more and better tests. In particular, if the user sets a CRS and the file contains a crs the projection is set based on what the user set but a warning is given that it doesn't agree with the file crs and the file CRS is given in the warning. It will now warn you in the situation above. Thanks for testing.

[shospital]

@rmendels that's great. I will check it out for the projection. I can also test other functions and get back to you.

[rmendels]

@shospital I have made a small change to xy_to_latlon() - it now uses the extract CRS if there no matter what the user gives. Will only use user CRS if none in file

[shospital]

@rmendels what if the metadata are incorrect? can user override it?

[rmendels]

@shospital Nope. If the CRS info from the data provider is incorrect, there are major problems which need to be corrected and which I do not want to touch. I am more concerned that the user gives a wrong CRS. I am going to stay with the data provider is correct.