From 3bf0f0eece1013e998f0c2ac1b05e50e9974e67e Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 6 Jun 2024 20:44:22 -0400 Subject: [PATCH] Fix all example errors --- examples/gmt/antioquia.jl | 3 ++- examples/is_it_a_plane.jl | 11 +++++++---- examples/multiple_crs.jl | 14 +++++++------- src/projection.jl | 10 ++++++---- 4 files changed, 22 insertions(+), 16 deletions(-) diff --git a/examples/gmt/antioquia.jl b/examples/gmt/antioquia.jl index 9844e274..ee7acaf3 100644 --- a/examples/gmt/antioquia.jl +++ b/examples/gmt/antioquia.jl @@ -29,6 +29,7 @@ f # using GeometryOps.jl to calculate the centroid. import GeometryOps as GO cx, cy = GO.centroid(antioquia_geoms) +# a.dest = "+proj=ortho +lon_0=$cx +lat_0=$cy" f # That looks a lot more like what the GMT example does! @@ -36,5 +37,5 @@ f # # make cover image #jl mkpath("covers") #hide -save("covers/$(splitext(basename(@__FILE__))[1]).png", fig) #hide +save("covers/$(splitext(basename(@__FILE__))[1]).png", f) #hide nothing #hide \ No newline at end of file diff --git a/examples/is_it_a_plane.jl b/examples/is_it_a_plane.jl index 04fd93ae..7320fe72 100644 --- a/examples/is_it_a_plane.jl +++ b/examples/is_it_a_plane.jl @@ -29,9 +29,9 @@ f, a, p = lines(reverse(Proj.geod_path(geod, reverse(jfk)..., reverse(sin)...)). # relative distances along the path, and altitudes of observation, as # a function of time. This is done by using the Animations.jl library. -times = [0, 0.5, 17.5, 18] +times = [0, 0.5, 30.5, 31] distances = [0, 0.05, 0.95, 1] -altitudes = [357860, 35786000/2, 35786000/2, 357860] +altitudes = [357860*2, 35786000/2, 35786000/2, 357860*2] distance_animation = Animation(times, distances, linear()) altitude_animation = Animation(times, altitudes, sineio()) @@ -47,11 +47,14 @@ satview_projection = lift(sl.value) do alt end ga = GeoAxis(fig[1, 1]; dest = satview_projection) meshimage!(ga, -180..180, -90..90, GeoMakie.earth(), shading = NoShading) +lines!(ga, GeoMakie.coastlines()) fig +# -record(fig, "plots/plane.mp4", LinRange(0, 1, 120)) do i - satview_projection[] = "+proj=geos +h=$(round(Int, at(altitude_animation, i*18))) +lon_0=$(Proj.geod_position_relative(inv_line, i)[2]) +lat_0=$(Proj.geod_position_relative(inv_line, i)[1])" +record(fig, "plane.mp4", LinRange(0, 1, 240)) do i + current_position = Proj.geod_position_relative(inv_line, i) + satview_projection[] = "+proj=geos +h=$(round(Int, at(altitude_animation, i*18))) +lon_0=$(current_position[2])" yield() end # ![](plane.mp4) diff --git a/examples/multiple_crs.jl b/examples/multiple_crs.jl index 192ff5bf..5c73e3e9 100644 --- a/examples/multiple_crs.jl +++ b/examples/multiple_crs.jl @@ -5,17 +5,17 @@ using CairoMakie, GeoMakie using Rasters, RasterDataSources, ArchGDAL CairoMakie.activate!(px_per_unit = 4) # hide -worldclim_temp = Raster(WorldClim{Climate}, :tmax; month = 1) -# Let's simulate a new CRS, assuming this was an image taken from a satellite: -new_worldclim = Rasters.warp( - worldclim_temp, +ras = Raster(EarthEnv{HabitatHeterogeneity}, :homogeneity) +# Let's simulate a new CRS, assuming this was an image taken from a geostationary satellite, hovering above 72° E: +projected_ras = Rasters.warp( + ras, Dict( - "s_srs" => convert(GeoFormatTypes.ProjString, Rasters.crs(worldclim_temp)).val, # source CRS + "s_srs" => convert(GeoFormatTypes.ProjString, Rasters.crs(ras)).val, # source CRS "t_srs" => "+proj=geos +h=3578600 +lon_0=72" # the CRS to which this should be transformed ) ) # This is what the raster would look like, if it were taken directly from a satellite image: -heatmap(new_worldclim; axis = (; aspect = DataAspect())) +heatmap(projected_ras; axis = (; aspect = DataAspect())) # Now, we can create a GeoAxis with coastlines in the equal earth projection: fig = Figure() ga = GeoAxis(fig[1, 1]) @@ -24,7 +24,7 @@ fig # The coastlines function returns points in the (lon, lat) coordinate reference system. # We will now plot our image, from the geostationary coordinate system: -surface!(ga, new_worldclim; shading = NoShading, source = convert(GeoFormatTypes.ProjString, Rasters.crs(new_worldclim)).val) +surface!(ga, projected_ras; shading = NoShading, source = Rasters.crs(projected_ras)) fig # Success! You can clearly see how the raster was adapted here. # diff --git a/src/projection.jl b/src/projection.jl index 7cc6a85f..f0b43f1d 100644 --- a/src/projection.jl +++ b/src/projection.jl @@ -140,10 +140,11 @@ end # GeoFormatTypes integration # ######################################## +const _GFTCRS = Union{GeoFormatTypes.CoordinateReferenceSystemFormat, GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS}} # Define methods for GeoFormatTypes CRS objects and all possible combinations thereof. -create_transform(dest::GeoFormatTypes.CoordinateReferenceSystemFormat, source::GeoFormatTypes.CoordinateReferenceSystemFormat) = create_transform(gft2str(dest), gft2str(source)) -create_transform(dest::String, source::GeoFormatTypes.CoordinateReferenceSystemFormat) = create_transform(dest, gft2str(source)) -create_transform(dest::GeoFormatTypes.CoordinateReferenceSystemFormat, source::String) = create_transform(gft2str(dest), source) +create_transform(dest::_GFTCRS, source::_GFTCRS) = create_transform(gft2str(dest), gft2str(source)) +create_transform(dest::String, source::_GFTCRS) = create_transform(dest, gft2str(source)) +create_transform(dest::_GFTCRS, source::String) = create_transform(gft2str(dest), source) """ gft2str(crs)::String @@ -152,4 +153,5 @@ Return a PROJ-compatible string from a GeoFormatTypes CRS object. """ function gft2str end gft2str(crs::GeoFormatTypes.EPSG{1}) = String("EPSG:$(GeoFormatTypes.val(crs))") -GeoFormatTypes2str(crs::GeoFormatTypes.CoordinateReferenceSystemFormat) = string(GeoFormatTypes.val(crs)) +gft2str(crs::GeoFormatTypes.CoordinateReferenceSystemFormat) = string(GeoFormatTypes.val(crs)) +gft2str(crs::GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS}) = string(GeoFormatTypes.val(crs))