Map Subsection USA States using matlab

I had the regional PM2.5 map, but I wanted to have it by states. So, came up with this code which can generate state level maps. If you use it, please let me know. Current example uses the solar radiation from GCM.

Highlight one, highlight all!

There is also R version available, at:

% Create the northeast mask
% Created by NM
% August 4, 2014
%
% http://www.mathworks.com/help/map/_f7-10852.html#f7-9156
% Use shaperead to get the patch data for the boundary:
close all
% http://www.mathworks.com/help/map/understanding-raster-geodata.html#f20-12849
pcs = {'New York', 'New Jersey','Massachusetts', 'Vermont','Connecticut' ,...
     'North Carolina','Ohio','Maine', 'New Hampshire', 'Rhode Island', ...
    'Maryland', 'Delaware', 'Virginia','West Virginia','District of Columbia' };
% pcs = {     'Pennsylvania'  };

cUS = shaperead('usastatelo.shp',...
    'UseGeoCoords', true,...
    'Selector',{@(name)any(strcmpi(name,pcs),2), 'Name'});

inLat = [cUS.Lat];
inLon = [cUS.Lon];
%%
% nys = shaperead('usastatehi.shp',...
%     'UseGeoCoords', true,...
%     'Selector', {@(name)strcmpi('New York',name), 'Name'});
% inLat = nys.Lat;
% inLon = nys.Lon;

% Set the grid density to be 40 cells per degree, and use vec2mtx to rasterize the boundary and generate a referencing vector for it:
gridDensity = 100;
[inGrid, inRefVec] = vec2mtx(inLat, inLon, gridDensity);
whos

% set the limits
figure(1)
axesm eqdcyl
meshm(inGrid, inRefVec)
colormap jet(4)

[latlim, lonlim] = limitm(inGrid, inRefVec);
setm(gca, 'Flatlimit', latlim, 'FlonLimit', lonlim)
tightmap

% To fill (recode) the interior of Indiana, you need a seed point
% (which must be identified by row and column) and
% a seed value (to be allocated to all cells within the polygon).
% Select the middle row and column of the grid and choose an index value of 3
% to identify the territory when calling encodem to generate a new grid:
inPt = round([size(inGrid)/2, 3]);
inGrid3 = encodem(inGrid, inPt,1);
% redraw the map using the filled grid:
inGrid3= 1-inGrid3; % need for multiple states

meshm(inGrid3, inRefVec)


%% Now try the Data
figure(2);
udGrid = flipud(inGrid3);
imagesc(udGrid)

[m,n] = size(udGrid);


% latlim=[36 46];
%     lonlim=[-82 -69];
     [Mlat, Mlon] = ndgrid(linspace(min(inLat),max(inLat),m), linspace(min(inLon),max(inLon),n));   %Changed from 35 46, -85 -69
%     newCF = griddata(latmatf,lonmatf,flipud(data),Dlat, Dlon);
 
    figure(3)
    ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
%     geoshow(ax, states, 'FaceColor', [1 1 1]);
    surfacem(Mlat,Mlon,(inGrid3));
 
 
    %% Now try the data
  load('C:\Users\CREST\Copy\CCNY-14\downscale-now\Regression-2006-1.mat')

    %     newCF = griddata(latmatf,lonmatf,flipud(data),Dlat, Dlon);
figure(4)
  ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Dlat,Dlon,isimipAve);
 
    %%
        newCF = griddata(Dlat,Dlon,isimipAve,Mlat, Mlon);
     
        inGrid3c = inGrid3;
        inGrid3c(inGrid3c<=0)= NaN;
        inGrid3c(inGrid3c>0) = 1;
        plotCF = inGrid3c.*newCF;
       
        figure(5)
     ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Mlat,Mlon,plotCF);
 
 
%% For Penn because When added Penn, it "inverts". So adding two with nanmean
pcs = {   'New York', 'Pennsylvania',  'Massachusetts', 'Vermont','Connecticut' ,...
     'North Carolina','Ohio','Maine', 'New Hampshire', 'Rhode Island', ...
    'Maryland',  'Virginia','West Virginia','District of Columbia' };

cUS = shaperead('usastatelo.shp',...
    'UseGeoCoords', true,...
    'Selector',{@(name)any(strcmpi(name,pcs),2), 'Name'});

inLat = [cUS.Lat];
inLon = [cUS.Lon];

% nys = shaperead('usastatehi.shp',...
%     'UseGeoCoords', true,...
%     'Selector', {@(name)strcmpi('New York',name), 'Name'});
% inLat = nys.Lat;
% inLon = nys.Lon;

% Set the grid density to be 40 cells per degree, and use vec2mtx to rasterize the boundary and generate a referencing vector for it:
gridDensity = 100;
[inGrid, inRefVec] = vec2mtx(inLat, inLon, gridDensity);
whos

% set the limits
% figure(1)
% axesm eqdcyl
% meshm(inGrid, inRefVec)
% colormap jet(4)

% [latlim, lonlim] = limitm(inGrid, inRefVec);
% setm(gca, 'Flatlimit', latlim, 'FlonLimit', lonlim)
% tightmap

% To fill (recode) the interior of Indiana, you need a seed point
% (which must be identified by row and column) and
% a seed value (to be allocated to all cells within the polygon).
% Select the middle row and column of the grid and choose an index value of 3
% to identify the territory when calling encodem to generate a new grid:
inPt = round([size(inGrid)/2, 3]);
inGrid3 = encodem(inGrid, inPt,1);
% redraw the map using the filled grid:
inGrid3=   inGrid3; % need for multiple states

meshm(inGrid3, inRefVec)
%

  inGrid3c = inGrid3;
        inGrid3c(inGrid3c<=0)= NaN;
        inGrid3c(inGrid3c>0) = 1;
        plotCF2 = inGrid3c.*newCF;
  figure(6)
     ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Mlat,Mlon,plotCF2);
 
    plotCF3d(:,:,1) = plotCF;
    plotCF3d(:,:,2) = plotCF2;
 
    iplotCF = nanmean(plotCF3d,3);
 
      figure(7)
     ax = worldmap(latlim, lonlim);
    states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
    geoshow(ax, states, 'FaceColor', [1 1 0.99]);
 
    surfacem(Mlat,Mlon,iplotCF);