Question about building a mesh with a specified CRS

[Houaga Isidore]

Dear Colleagues,

I trust you are fine. I’m having a challenge in building a mesh with a desired CRS using the script below.

 

I believe that the issue is because the map coordinates are in degrees while my desired crs (crs_tz) is not.

 

Please, how can I transform the Tanzania map to match my desired CRS. Thanks. 

 

Best Regards,

Isidore

 

#--------------------------------Building a new mesh --------------------

# Step1: Get map of Tanzania

mapTZA <- getData('GADM', country = "TZA", level= 1) # Get map of Tanzania

# Step 2: Define the desired CRS (EPSG:21035 - Arc 1960 / UTM zone 35S)

crs_tz <- fm_crs("epsg:21035")  # UTM 35S for Tanzania

fm_length_unit(crs_tz) <- "km"  # Set units to kilometers

 

# Step 3: Check and transform mapTZA to the desired CRS if necessary

if (!identical(proj4string(mapTZA), proj4string(crs_tz))) {

  # Transform mapTZA to EPSG:21035

  mapTZA <- spTransform(mapTZA, crs_tz)

}

 

# Step 4: Use gUnaryUnion to merge the regions into a single boundary for Tanzania

TanzaniaBorder <- gUnaryUnion(mapTZA, rep(1, nrow(mapTZA)))

 

# Check and print the CRS of TanzaniaBorder to confirm

print(proj4string(TanzaniaBorder))

 

# Extract the border from the map

TanzaniaBorder <- gUnaryUnion(mapTZA, rep(1, nrow(mapTZA)))

# Formatting boundary for r-INLA

TanzaniaBorder <- inla.sp2segment(TanzaniaBorder)

TanzaniaBorder <- fm_as_segm(TanzaniaBorder)

plot(TanzaniaBorder)

bnd <- fm_nonconvex_hull_inla(locations, 0.5)

 

gps_data <- st_as_sf(data1, coords = c("long", "lat"), crs = 4326)  # CRS: WGS 84 (EPSG: 4326)

transformed_data <- fm_transform(gps_data, crs = crs_tz)

locations <- st_coordinates(transformed_data)

 

mesh <- fm_mesh_2d_inla(boundary = list(bnd, TanzaniaBorder),

                        max.edge = c(12.8, 64), offset=5.2, cutoff = 25, crs = crs_tz)

# Plotting mesh

ggplot() +

geom_fm(data = mesh) +

geom_point(aes(locations[,1], locations[,2])) + labs(x="Longitude", y="Latitude")

 

 

Fig1. Resulting mesh from the script above.

 

 

Fig 2. The expected mesh (I want to produce a mesh like this). This mesh was based on CRS: WGS 84 (EPSG: 4326).

Best Regards,

Isidore Houaga, PhD

Research Scientist

Science Manager at AABNet (https://www.animalbreeding-africa.org)

Secretary General AfricaBP (https://africanbiogenome.org/)

Royal Society Newton International Fellow (Alumni)

Highlander Lab  & Centre for Tropical Livestock Genetics and Health (CTLGH)

The Roslin Institute | University of Edinburgh Easter Bush Campus | Midlothian EH25 9RG | United Kingdom

 E-mail: isidore...@roslin.ed.ac.uk / iho...@gmail.com 

www.ctlgh.org  www.ed.ac.uk/roslin http://animalbreeding-africa.org/
                         

 The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

[Elias T. Krainski]

Hi,

I was not able to reproduce your code, and could not open the two figures. I have a code to make a mesh around a country (using a map retrieved with the rnaturalearth package). However, I didn't succeed in fetching a map of Tanzania. The attached code takes a map from nbs.go.tz to do it.

regards,

Elias

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAC1Lj_BewwWQ%2BU3FSTki2YTgxLbEi4QS4BiC43YY_V%3D5bgnfcg%40mail.gmail.com.