How to plot a contour map in UTM projection?
9 views (last 30 days)
Show older comments
I am trying to plot data over a regular grid(GRID_VEL50 corresponding to Northings and Eastings in meters-UTM map projection) in the form of filled contours. I want the corresponding map overlayed along with it. Following are the variables used in the code: GRID_VEL50 (dim:156X150),GRID_X (Eastings,dim:150) and GRID_Y (Northings,dim:156). Here is the script, I tried for plotting the same:
axesm('tranmerc','MapLatLimit',[7 37.5],'MapLonLimit',[66.667 98.6666717529297]);
framem; gridm; axis off; tightmap;
x11 =245895.96; %meter, the North-West corner(66.667 deg in long)
y11 =774370.68; %meter, the North-West corner(7 deg in lat)
dx1 = 25000; %in meter
dy1 =-25000; %in meter
R=makerefmat(x11, y11, dx1, dy1);
contourfm(GRID_VEL50,R,'r','LineWidth',0.0001)
coast = load('coast');
geoshow(coast.lat, coast.long, 'Color', 'black')
contourcbar
On running this script, it gives an error, stating :"Error using constructGeoRasterReference (line 45)." The referencing array comes out to be a 3X2 array as mentioned in help documentation.
R =
0 -25000
25000 0
220895.96 799370.68
I am unable to debug this. Is there any other method to plot contours in a specific map projection?
4 Comments
Walter Roberson
on 10 Sep 2015
You indicate lat 7 long 66.667 as your North-West corner, but your axesm latitudes are from 7 to 37.5 which implies that lat 7 is to be a South corner rather than a North corner. Please re-check which points are the North or South extremes and please check the sign of your dx1 and dy1 .
When I calculate with the dx1 and dy1 values you give, I do not as yet see any reason why the latitude would exceed +90 to -90, but I do see that it would be outside your map axes which suggests you have incorrect bounds or incorrect signs.
Answers (1)
Walter Roberson
on 11 Sep 2015
In the documentation it is not clear to me how makerefmat distinguishes between the x/y case and the long/lat case. I speculate that if you pass integral numbers that it might interpret them as x/y, so try with
R=makerefmat(round(x11), round(y11), round(dx1), round(dy1));
See Also
Categories
Find more on Map Display in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!