Overlaying a grid of colored squares onto a map in GMT

[yalonso]

Hello! I am new to Julia and GMT and am looking to see if it’s possible to overlay a grid (in which i define how large each box should be and the color it should be) on top of a plot of a certain region of the US.

function create_grid(region, grid_size)
    lon_min, lon_max, lat_min, lat_max = region
    lon_step = (lon_max - lon_min) / grid_size
    lat_step = (lat_max - lat_min) / grid_size

    grid = [((lon_min + i*lon_step, lat_min + j*lat_step), 
        (lon_min + (i+1)*lon_step, lat_min + (j+1)*lat_step)) 
            for i in 0:(grid_size-1), j in 0:(grid_size-1)]
    return grid
end

I have written this function to define how large the “Squares” of the grid should be (it just contains the bottom left and top right points of each square). I have values for each square that would help me define what the color of the square should be. I thought that using grdimage (grdimage · GMT) could be helpful, but I’m not entirely sure how to create a .grd file that contains the boundary information of each square. Is it possible to define the grid by the vertices of each square? Thank you so much. I’ve attached a rough sketch of what I’d like for the plot to look like!

Update: we are now trying to use a method similar to what is described here and this is our code to do so:

function make_plot_with_grid()
    # Define the latitude and longitude of Hayward, CA
    hayward_lat = 37.6688
    hayward_lon = -122.0808
    
    # Define the zoom level
    zoom_level = 0.5

    # Calculate the region coordinates
    my_region = [hayward_lon - zoom_level, hayward_lon + zoom_level, hayward_lat - zoom_level, hayward_lat + zoom_level]
    lon_min, lon_max, lat_min, lat_max = my_region

    grid_size = 20
    lon_step = (lon_max - lon_min) / grid_size
    lat_step = (lat_max - lat_min) / grid_size
    
    fig = coast(region=my_region, proj=:Mercator, frame=:g, area=1000, land=:gray, water=:white)
    colors = zeros(grid_size, grid_size)
    for i in 1:grid_size
        for j in 1:grid_size
            sub_lon_min = lon_min + i * lon_step
            sub_lat_min = lat_min + j * lat_step
            sub_lon_max = lon_min + (i + 1) * lon_step
            sub_lat_max = lat_min + (j + 1) * lat_step
            cell = (sub_lon_min, sub_lat_min), (sub_lon_max, sub_lat_max)

# prob gets a value that we base the color of the cell off of
            prob = sample_and_average(cell, myMDP, 10)
            color = probability_to_color(prob)
            colors[i, j] = color
        end
    end

    coast!(region=my_region, proj=:Mercator, figsize=50) # make the map
    topo = makecpt(color=:rainbow, range=(0,225,20), continuous=true) #make colorbar

    grdimage!(colors, coast=true, color=topo, limits=my_region)#frame=(annot=30, grid=30, ticks=20, title="Our test title"), par=(FONT_ANNOT_PRIMARY=18,)) # plot my gridarray
    #colorbar!(pos=(anchor=:BC,length=(32.5,0.7), horizontal=true, offset=(0,1.0)),color=topo, frame=(ylabel="W/m^2",),par=(FONT_ANNOT_PRIMARY=18,), show=true,savefig="/home/yasminealonso/Desktop/gmt1.png") # plot colorbar & show and save figure
    # Display the plot
    show(fig)
end

but we are getting this error:

grdimage [WARNING]: Your grid y's or latitudes appear to be outside the map region and will be skipped.
grdimage [WARNING]: No grid or image inside plot domain
nothing

Thank you so much!

[joa-quim]

Don’t have much time right now to read this with attention but I think you want to use pcolor followed by a coast! command to add the coastlines. Please try and come back if it doesn’t solve it for your case.

[yalonso]

Thank you – I think this definitely could work! For some reason though, the area using pcolor only spans a very small area over the region that I define in coast. I’ve printed out the coordinates and it does seem like the coordinates I’ve put into the X and Y vectors should completely span the entire coast region. Have you seen this before? Thank you so much for your helpful and timely response!

[joa-quim]

Hmm, not sure how you got that but here are two alternatives (with a dumb grid):

x = -122.6:0.1:-121.8; y = 37.1:0.1:38.0;

julia> [length(y) length(x)]
1×2 Matrix{Int64}:
 10  9

# Note how the matrix size is 1 shorter than vector lengths
mat = rand(Float32, 9,8);

# You have to plot the tiles first otherwise they would hide coast lines
pcolor(x, y, mat, proj="merc")
coast!(limits=[-122.6,-121.8, 37.1, 38.0], show=true)

But you probably seek something more on the line of

# Use the same x,y and mat as above
G = mat2grid(mat, x=x, y=y, reg=1, proj="lonlat")

grdimage(G, proj=:merc, coast=true, show=true)

[yalonso]

Thank you so much – this is so helpful. I just have one last question – do you know if it is possible to mask out the grid from places in the water? I’d like to have the colored squares only over land and have the water area just be grey (not gridded up). Thanks for all your advice!

[joa-quim]

grdlandmask let’s you compute a grid mask that you can multiply or grid with. mask let’s you do image clipping. But see also sealand and this gallery example