Alternative to mapprofile function for MapCellsReference
10 views (last 30 days)
Show older comments
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!
0 Comments
Answers (0)
See Also
Categories
Find more on Geographic Coordinate Reference Systems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!