# Projecting and interpolating data to a regular grid

**URL:** <https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162>\
**Category:** Geo\
**Tags:** interpolations, gmt, geodesy\
**Created:** [February 27, 2021, 10:38pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162 "2021-02-27T22:38:41Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![johnbb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnbb/32/34233_2.png) [@johnbb](https://discourse.julialang.org/u/johnbb)\
**Post date:** [February 27, 2021, 10:38pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/1 "2021-02-27T22:38:41Z")

</div>

I have (multiple) data on a regular latlon grid which I want to project with a Lambert conformal conic projection and interpolate to a regular grid. Is using the Proj4 and Interpolations packages an appropriate option or is there a more tailored alternative?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 28, 2021, 6:31pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/2 "2021-02-28T18:31:28Z")

</div>

Would it be possible to provide some more details, some reference or some images illustrating your problem?

I am not an expert in Geodesy, just an interested user, and in the mock-up example below produced with [Proj4.jl](https://github.com/JuliaGeo/Proj4.jl), I have tried to visualize your problem, probably in a wrong way, but hope it helps getting more feedback.

The left image is a regular Lat-Lon grid in WGS84 and the right image its transformation to North America Lambert conformal conic (LCC) X-Y coordinates.

It does not make immediate sense to define a regular grid on this output. What would be the application?

_PS: from @joa-quim’s GMT grdproject reference provided, one may define a regular grid that fits inside the projected LCC irregular grid, then do an inverse projection and perform bicubic spline interpolation (ex: with Dierckx.jl Spline2D) using the available spatial data (ex: temperature T(X,Y)) in the original grid._

 ![Lambert_interpolation](https://global.discourse-cdn.com/julialang/original/3X/9/0/903a3eabd595c8dc73c497a47997243e5408097b.png)

---

<div class="post-metadata">

**Author:** ![johnbb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnbb/32/34233_2.png) [@johnbb](https://discourse.julialang.org/u/johnbb)\
**Post date:** [February 28, 2021, 7:10pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/3 "2021-02-28T19:10:11Z")

</div>

Yes, your example illustrates the first part well. The next step is to interpolate the data on the Lambert grid to a grid expanded by the vectors, say,  
x = 1e6:1e4:2e6  
y = 3e5:1e4:2e6  
The aim is to match two datasets with different projections.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [February 28, 2021, 7:17pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/4 "2021-02-28T19:17:15Z")

</div>

If you have a grid in geographical coordinates you can reproject that grid into Lcc and get regular cartesian grid in _Lcc meters_. Either GMT’s `gdproject` or GDAL `gdalwarp` let you do that.

---

<div class="post-metadata">

**Author:** ![johnbb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnbb/32/34233_2.png) [@johnbb](https://discourse.julialang.org/u/johnbb)\
**Post date:** [February 28, 2021, 7:51pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/5 "2021-02-28T19:51:42Z")

</div>

Thanks. I will need to look further into these. Are there any Julia examples of using `GMT.grdproject`?

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [February 28, 2021, 8:13pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/6 "2021-02-28T20:13:38Z")

</div>

For GMT you would need to create a GMTgrid type. I’m sorry but cannot give a pointer to the docs because something screw in the documentation builder and the GMT docs are giving a 404 right now. Hope it recovers soon.  
But the GMT’s `mat2grid` function offers lots of possibilities to create them. It has an online help.

The, assuming **G** holds the geog grid, the command would be something like (need to complete the parameters for projection using the PROJ4 syntax)

Glamb = grdproject(G, proj=“+proj=lcc+lat\_1=…+lat\_2=…+lat\_0=…”) ;

and eventually even easier if you already have the geog grid as a file on disk. Then the comand would be (assuming that GMT will guess the file format).

Glamg = grdproject(“filename”, proj=“+proj=lcc+lat\_1=…+lat\_2=…+lat\_0=…”) ;

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 28, 2021, 10:43pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/7 "2021-02-28T22:43:28Z")

</div>

@johnbb, if the regular grid is larger than the projected input data, extrapolation will be also required. There will be some geometrical/scale distortion. See also the postscript in my updated post above for additional comments.

---

<div class="post-metadata">

**Author:** ![johnbb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnbb/32/34233_2.png) [@johnbb](https://discourse.julialang.org/u/johnbb)\
**Post date:** [February 28, 2021, 11:04pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/8 "2021-02-28T23:04:13Z")

</div>

> [@rafael.guerra](#):
>
> PS: from @joa-quim’s GMT grdproject reference provided, one may define a regular grid that fits inside the projected LCC irregular grid, then do an inverse projection and perform bicubic spline interpolation (ex: with Dierckx.jl Spline2D) using the available spatial data (ex: temperature T(X,Y)) in the original grid.

The inverse approach is appealing as I have many datasets and the latlon grid is less dense, but covers completely the target area/grid. Proj4 and Interpolations may then be a good alternative.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [February 28, 2021, 11:04pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/9 "2021-02-28T23:04:44Z")

</div>

OK, a full GMT example

```julia
# Extract a sub-region from a global grid. Using a low resolution to keep example light
Ggeog = grdcut("@earth_relief_10m", region=(-18,15,32,56));

```

convert to Lambert

```julia
Glamb = grdproject(Ggeog,proj="+proj=lcc+lat_1=39+lat_2=50+lat_0=45+lon_0=0");

```

Show them. Rafael, you are right about the extrapolation but the trick is to not do it and instead fill those nodes with NaNs.

```julia
imshow(Ggeog, shade=true, cmap=:geo, name="geogs.png")
imshow(Glamb, shade=true, cmap=:geo, name="lcc.png")

```

 ![geogs](https://global.discourse-cdn.com/julialang/original/3X/6/a/6a1f7f369ac21f9537a19e0d2f48900e29f17416.png)  
 ![lcc](https://global.discourse-cdn.com/julialang/original/3X/f/2/f2f48dab880845236657b72d7fa95613cdd3c9bf.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 28, 2021, 11:11pm UTC](https://discourse.julialang.org/t/projecting-and-interpolating-data-to-a-regular-grid/56162/10 "2021-02-28T23:11:31Z")

</div>

@joa-quim brilliant plot. Love it.