How can I georeference this grid/matrix?

4 次查看（过去 30 天）
Harith Kubaisy
Harith Kubaisy 2020-7-26
评论： Maskus Kit 2021-2-6
Hello,
I have a script which plots kernel density values for earthquakes. However, I need to visualize the final results of the p array on a proper projected geograhpic map. How can I incorporate the results on a georeferenced MATLAB map? If that is not possible, what is the best way to export the data to ArcGIS platform?
I tried saving the results as a .txt file (ASCII) and performed a conversion to Raster by ArcMap 10.4.1 but I could not get meaningful results.
The script can be seen below, and the longitude/latitude files are attached.
clc,clear
load long5
load lat5
long = (long5);
lat = (lat5);
step = 0.009;
x = min(long):step:max(long);
y = min(lat):step:max(lat);
[X,Y]= ndgrid(x,y); % creates grid
R = 0.05; % circle radius
plot(long,lat,'.r') % plot data
P = X*0; % no. of points inside each circle
T = 2019-1916+1;
sig = 2;
scale = (pi/180*earthRadius/1000)^2;
hold on
for i= 1:length(x)
for j=1:length(y)
D2 = (long-x(i)).^2 + (lat-y(j)).^2;
ix = D2 < R^2; % find smaller distances
D(ix) = scale*D2(ix); % write smaller distances
P(i,j) = sum(exp(-D(ix)/2/sig^2))/T; %write ID
N(i,j) = sum(ix); % no. of points inside (i,j) circle
% plot(long(ix),lat(ix),'ob') %plot data inside circle (if exists)
if sum(ix) % if there are points inside
viscircles([x(i),y(j)], R,'edgeColor','g'); % plot circle
else
% viscircles([x(i),y(j)], R,'edgecolor','r'); % plot circle
end
end
end
p = P/(pi*6.^2);
G = rot90(flip(p),3);
imagesc(x,y,G);
colorbar
hold off
axis equal
  3 个评论
显示 1更早的评论
Harith Kubaisy
Harith Kubaisy 2021-2-5
Hello. Yes the following work-around did it for me. Sorry for not updating the post.
In a seprate script:
V1 = X(:);
V2 = Y(:);
PV = pMo(:);
lon = (V1);
lat = (V2);
weights = (PV);
hold off
figure
geodensityplot(lat,lon,weights,'FaceColor','interp');
hold on
geobasemap grayterrain
colormap(hsv)
colorbar

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

jonas
jonas 2020-7-26
With the mapping toolbox you can just replace your plot functions with the corresponding map axes function.
load coastlines
ax = worldmap('italy');
plotm(coastlat,coastlon)
plotm(long,lat,'.r') % plot data
h = pcolorm(x,y,G);
should give you something like this
  1 个评论
Harith Kubaisy
Harith Kubaisy 2020-7-27
Thank you for the quick response, I tried the approach but I did not get the plot over my desired area.
See attached script and result map.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Geographic Plots 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by