Alternative to mapprofile function for MapCellsReference

10 views (last 30 days)
I am trying to use the mappin toolbox to generate an elevation profile along a prescribed path. I have elevation data as a GeoTIFF file. I have imported elevation data from a .tif file with the function readgeoraster, and the path from a .gpx file with the function gpxread, therefore I have a matrix A with elevation data (in meters), a geopoint vector path with info of my path (in latitude and longitude), and a geographic cells reference structure R with info about the coordinates to which A is referred. Then I use the function mapprofile to retrieve elevation elev along the path and also distance dist between the interpolated points. Here is my code:
close all; clear; clc;
[A,R] = readgeoraster("output_COP30.tif","OutputType","double");
% .tif file downloaded from https://portal.opentopography.org/datasets
% copernicus 30 m dataset follows the standard WGS 84 (i.e., it is in "degrees")
routeData = gpxread("route.gpx");
% .gpx file downloaded from https://map.project-osrm.org/
[elev,dist,lat,lon] = mapprofile(A,R,routeData.Latitude,routeData.Longitude);
However, when I import a GeoTIFF file that has been made as a projection (in my case, " UTM zone 32N") I obtain a map cells reference structure as R, and now I cannot use the function mapprofile anymore. Here is the alternative code:
close all; clear; clc;
[A,R] = readgeoraster("w45010_s10.tif","OutputType","double");
% .tif file downloaded from https://tinitaly.pi.ingv.it/data_1.1/w45010_s10/w45010_s10.zip
% TinItaly dataset is projected, via WGS 84 / UTM zone 32N (i.e., it is planar)
routeData = gpxread("route.gpx");
% .gpx file downloaded from https://map.project-osrm.org/
[elev,dist,lat,lon] = mapprofile(A,R,routeData.Latitude,routeData.Longitude);
I need to use the second dataset because it has a high resolution (10 m instead of 30 m of the first one), which better suit my application.
Is there a way to do the same operations that mapprofile performs but on a "projected" reference system?
Alternatively, is there a way to "unproject" / interpolate / translate my "projected" data in order to obtain a geographic cells reference and still use mapprofile?
I also tried to generate a new geographic cells reference by using georefcells and then use it in mapprofile, as follows:
R_info = georasterinfo("w45010_s10.tif");
[latlim,lonlim] = projinv(R_info.CoordinateReferenceSystem,R.XWorldLimits,R.YWorldLimits);
R_geo = georefcells(latlim,lonlim,R.RasterSize,'ColumnsStartFrom',R.ColumnsStartFrom);
[elev_1,dist_1,lat_1,lon_1] = mapprofile(A,R_geo,routeData.Latitude,routeData.Longitude);
but it does not seem to retrieve correct elevation and it is completely wrong about distance (because it misses the GeographicCRS property). Moreover, I guess this is not a proper "unprojection" because only latitude and longitude limits (i.e., vertices of the square used) are correctly translated, but the internal values of the matrix A will not refer to the "unprojected" points, since the grid will not be made in the same way.
For simplicity, I added the two .tif files and the .gpx file in the shared Drive folder: https://drive.google.com/drive/folders/1jj1KIe7RJQTf3XzRow2vBpUOkyY-8TbK?usp=drive_link
Thanks for any suggestions!

Answers (0)

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!