How can I georeference this grid/matrix?

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 

[Maskus Kit]

Hey, did you finally figure this out?

[Harith Kubaisy]

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

[Maskus Kit]

Thanks a bunch. Really helpful!

Answer 1

[jonas]

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

[Harith Kubaisy]

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.