alternative projections

[smormede]

Hello. I was wondering if we could use our own projection rather than UTM (to match with other data we want to apply to the results) and how difficult that would be. Is it embedded through the code or only a few places where we’d have to change it?

I see rgdal has been mentioned in other threads, and this is what I have used in R, with different projections.

thanks

Sophie

[James-Thorson]

Do you intend to use a projection for which there is a function in R to convert to/from Latitude/Longitude to you projected coordinates? Or do you instead want to simply input the projected coordinates for your observations and the extrapolation grid, without specifying a function in R to do the conversion?

[smormede]

Hello. I can do either, it depends if the functions inside delta-glmm or VAST need to know the transform or not. I would most likely transform in R using rgdal, the point being to have the same projection as other pieces of analysis I want to use in conjunction with this one. Here is an example with the projections I have used as an example:

library(rgdal) dstart<-data.frame(lon=dat2$long_sj, lat=dat2$lat_sj) # that's the object coordinates(dstart) <- c("lon", "lat") proj4string(dstart) <- CRS("+init=epsg:4326") # that's the lat long projection CRS.new <- CRS("+init=epsg:2193") # that's the eastings and northings projection dstart.t <- spTransform(dstart, CRS.new) # here's where you transform

dstart.t@coords[,"lon"] # longitude in the projection

[James-Thorson]

Hmm. This looks like it should be useful, and it should be feasible to add an option for using this projection method instead of the PBSmapping package for UTM projections. However, I'd like to learn more about rgdal and the projections its using before committing to this path.

• Do you have a recommended reading to explain what projection is meant by epsg:4326 and epsg:2193?
• For example, why is the lat long projection not using +proj=longlat?\

I apologize, but reading ?spTransform is a bit dense!

[smormede]

Hi Jim

the EPSG projections are pre-defined and have the appropriate formula attached inside rgdal. These are defined internationally, see for example http://spatialreference.org/ref/epsg/

EPSG 4326 corresponds to WGS 84, the lat and long coordinates we use. So that line just says that the data I am giving is in lats and longs. https://en.wikipedia.org/wiki/World_Geodetic_System

EPSG 2193 is the equal area projection that gets used often around New Zealand and was used for other parts of the project I want to link with. If you click on "Proj4" in the link below you will get the exact projection definition, copied below the link. http://spatialreference.org/ref/epsg/2193/ +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

The value of this second projection would be the one that people a most likely to want to change for their specific area of interest. It can be given as an epsg number or self defined.

I hope this makes sense.

Sophie

[James-Thorson]

Sophie,

Thanks for helping me learn a bit more about spatial data standards, this looks very promising!

Am I correct in assuming that the CRS for the transform will always involve a character (e.g., "+init=epsg:2193")? If yes, I might overload the meaning of the zone input to my existing functions, such that if its integer it uses UTM and if its character it uses this new method.

Jim

[smormede]

Hi Jim. I believe so. That would be great. A couple more links that might be useful (if you haven't already found them): https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf

https://www.rdocumentation.org/packages/rgdal/versions/1.2-15/topics/CRS-class

[James-Thorson]

OK, I just pushed a commit to SpatialDeltaGLMM, which should be updated automatically if you re-install VAST from github. The new feature is accessed by changing code at a single location

Extrapolation_List = Prepare_Extrapolation_Data_Fn( Region=Region, strata.limits=strata.limits, zone="+init=epsg:2193" )

where zone previously was assumed to be an integer (UTM zone), and now if its a character it uses rgdal and interprets the character-input as the CRS. This feature is only built into the backend function Prepare_NZ_Extrapolation_Data_Fn so its not available for other regions, but would be easy to add.

I ran through the Chatham Rise hake example, and it all seemed to be working. However, I didn't do any testing on a branch (because SpatialDeltaGLMM is a dependency of VAST and I don't know how to use branching to test dependencies well). So please use cautiously and if you find a bug please use the minimal-example reporting https://github.com/nwfsc-assess/geostatistical_delta-GLMM/wiki/How-to-create-a-minimal-example. Also, if its working please close the issue.

[James-Thorson]

Just ran it on a clean terminal and caught a small bug (didn't export the new function to the namespace).

One other update: you'll need to manually update SpatialDeltaGLMM via devtools::install_github("nwfsc-assess/geostatistical_delta-GLMM")

[smormede]

Hi Jim,

The projection seems to work, even with the full formula, which is great. I get exactly the same locations in Spatial_List$loc_x if I calculate them externally.

However it looks like "flip_around_dateline" is still used and creating issues not only in the plotting but also the model itself somehow: the model results look flipped around 180 degrees, but also the cvs calculated are completely erroneous. If I set flip_around_dateline = F in Prepare_Extrapolation_Data_Fn then the model results are comparable with the model using UTM. I suggest either some instructions saying "set flip as F if you use rgdal", or coding it so it doesn't happen.

thanks

[smormede]

Update about the above for people interested in using their own projection. At the moment it's only available for the New Zealand Extrapolation_List. But you can still use it as follows, I haven't found other places to change it yet. Only the relevant changes are included below.

Region<-"New_Zealand"

Extrapolation_List = SpatialDeltaGLMM::Prepare_Extrapolation_Data_Fn(Region = Region,strata.limits = strata.limits,Data=temp.dat,zone="+proj=laea +lat_0=-90 +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0",flip_around_dateline=F)

here I am giving the Ross Sea projection as an example.

Get region-specific settings for plots

MapDetails_List = SpatialDeltaGLMM::MapDetails_Fn( Region = Region, NN_Extrap = Spatial_List$PolygonList$NN_Extrap, Extrapolation_List = Extrapolation_List )

need to hack into the map details for the x and y lims otherwise it gives the Chatham Rise region

MapDetails_List$Xlim<-c(150,210) MapDetails_List$Ylim<-c(-79,-60)

and then you can change the plotting ratio

source(make.filename("geodesicDistance.r",DIR$Functions)) #or whatever function you use to calculate distance temp<-DistanceLongLat(MapDetails_List[["Xlim"]][1],MapDetails_List[["Xlim"]][2],MapDetails_List[["Ylim"]][1],MapDetails_List[["Ylim"]][1])/DistanceLongLat(MapDetails_List[["Xlim"]][1],MapDetails_List[["Xlim"]][1],MapDetails_List[["Ylim"]][1],MapDetails_List[["Ylim"]][2]) MapDetails_List[["MapSizeRatio"]][1]<-MapDetails_List[["MapSizeRatio"]][2]/temp