New answers tagged raster
1
The expansion is done because of the coord_sf() parameter expand = TRUE.
Setting it to FALSE should remove the "buffer".
Here's the corresponding part of the documentation.
ggplot() +
geom_point(data= basemap.df, aes(x=x, y=y,col=rgb(layer.1/255, layer.2/255, layer.3/255))) +
scale_color_identity() +
geom_sf(data=state, color= &...
0
This is your decision.
I think there are two ways.
Calculate by interpolating the missing values.
The average is obtained from the corrected value.
The mean is calculated only when each cell has a value other than Nodata.
If you use only the QGIS raster calculator, first classify the values of each raster as 1 and nodata=0.
Adds all rasters.
You will be ...
0
show(stack[0])
show(stack[1])
will work fine.
If you choose tuple, you have to choose the channel to output. Multi-band output is possible only when the raster dataset is loaded.
The same is true for the creation of raster files.
dst.write(stack, 1)
will change
dst.write(stack[0], 1)
dst.write(stack[1], 2)
UPDATE
The cause of the error is that the width ...
2
Sanborn has an online image purchase portal at: https://online.sanborn.com/
Maybe they have coverage in your area with suffucient rexooution.
(Note that I am not affiliated with Sanborn and do not receive any benefit from this recommendation...except possibly SE reputation points!)
2
You are trying to read raster data with read_csv:
library(readr)
img1 <- read_csv("E:/Topic_sinkholes/Final data/stack layers/img1.hdr")
You've not told us where you got this file from but if it is geospatial raster data it is unlikely to be a CSV file.
Try:
library(raster)
img1 = raster("E:/Topic_sinkholes/Final data/stack layers/img1.hdr&...
1
you can simply crop one crop to another and end up with the intersection of all of them. try
r.all <- list(r1,r2,r3,r4,r5) # your 5 rasters
for (r in r.all) {
xmin <- c(NA,extent(r)[1])
xmax <- c(NA,extent(r)[2])
ymin <- c(NA,extent(r)[3])
ymax <- c(NA,extent(r)[4])
}
#bounding box of the intersection of all rasters
b.box &...
2
You can do this in raster by unioning extents. You don't need st_bbox. This way you get back an extent object ready for use in the raster package functions you might need:
> r1 = raster(matrix(1:12,3,4),xmn=.3,xmx=.5,ymn=.2,ymx=.8)
> r2 = raster(matrix(1:12,3,4),xmn=.4,xmx=.5,ymn=.3,ymx=.9)
> union(extent(r1), extent(r2))
class : Extent
xmin ...
0
You could call gdal_translate on the input raster and use the -projwin
argument to specify the bounds you require.
Since you know the position of the city x & y to calculate sides of your box is simple maths.
side = 30000 # 30km
h_side = side//2 # half a side (15km)
left = x - h_side
right = x + h_side
top = y + h_side
bottom = y - h_side
then simply ...
2
It's this bit, is wrong it assigns the crs doesn't apply it by transforming the data:
projection(landcover) <- CRS("+init=epsg:32735")
you're expected to transform the raster i.e.
new <- raster(healthfac_point_projected)
res(new) <- c(150000, 150000) ## change resolution to suit, this is metres (the width,height of each cell)
landcoverutm ...
1
You're writing a "MEM" raster, it's in MEMory. If you want to write to file, use a file based format like a GeoTIFF.
From the docstring in @Jose's rasterize_me function:
If you want to generate a GTiff on disk, set format to GTiff and fname_out to a sensible filename.
So do something like:
rasterize=rasterise_me(r"the/path/to/my/raster.tif&...
0
The raster package has a handy wrapper function for reading shapefiles. You should be able to do:
library(raster)
CEZ <- shapefile("./Shapefile/CEZ1.shp")
if that's where your shapefile really is. You might get a better error message.
0
For those of you who are interested in this topic, this document provides lots of information about different binning methods.
Other useful information are in the comments of this github page which contains very useful python implementations of L2 data binning.
2
Your input raster might not be filled with NODATA values but still, it might be defined to use a certain NODATA values when NODATA values are necessary.
You can check the NODATA value defined for your raster with the ST_BandNoDataValue(rast, band) function inside PostGIS or with gdalinfo outside the database.
If your raster have a different NODATA defined ...
1
The Previous Answers worked for QGIS 2.x.
If you are reading this after 2020, you will find that the Answers don't work. What Works with QGIS 3.x is the plugin called IBAMA Processing.
Install this plugin, and then you will get the `Create Footprint of Images' tool/algorithm in the Processing Toolbox.
1
Check out the vectorization settings within the ArcScan toolbar. There is a Maximum Line Width (MLW) setting that can be adjusted within a range of 1-100. Unfortunately, 100 is the maximum. FYI: a line width value is the product of the coordinate system measurement unit and the pixel size. If, at a MLW setting of 100, ArcScan doesn't vectorize a centerline ...
1
I faced this type of issue when importing and write ESA SNAP .dim images in R. Then after some research I tried to update the crs of the raster from here. It works for me.
try this for your case:
crs(raster_objects)<-CRS("+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs")
then export your image. you can find the proj4string it here.
2
Try "Paletted/Unique values" and color ramp "Greys" under symbology tab:
4
Starting with your raster stack s:
> s
class : RasterStack
dimensions : 15, 10, 150, 6 (nrow, ncol, ncell, nlayers)
I'll show how to fit and predict in various ways. I'm going to try to spell out every stage and use data structures that make it clear what's going on - some of these steps can be made quicker in various ways but I'm aiming for ...
0
Perhaps an alternative idea, instead of using the intersect function, a common way in QGIS/ArcGIS Pro on how to deal with such problems is to first calculate an intersection point vector at each of the road networks polygon boundary.
You could then use some ready made tools such as discussed here (SAGA - Split lines at points), to calculate/split the length ...
2
You probably need to enable the CUTLINE_ALL_TOUCHED warp option:
CUTLINE_ALL_TOUCHED: This defaults to FALSE, but may be set to TRUE to enable ALL_TOUCHEd mode when rasterizing cutline polygons. This is useful to ensure that that all pixels overlapping the cutline polygon will be selected, not just those whose center point falls within the polygon.
4
You can rasterize the map using Project->Import/Export -> Export Map to Image. You can check "Append georeference information" and no need to georeference. You can select from several raster formats (tiff, jpg, png, etc). Dont't forget to copy the world file (jgw, pgw, etc) to keep georeference.
0
This is not a full answer but may help someone with the same issue.
I suceeded in testing my raster file on a different installation and the file loads successfully.
Problem appears to be something wrong with my installations of 3.10 and 3.12.
As I have a working version 3.4, I will not risk fixing the problem at the moment. The update to OSTN15NTv2 did not ...
1
You can create a numpy Array using random choice, convert to a ascii raster then add to your dtm using raster calculator:
import numpy as np
proportion_1 = 0.1 #Adjust
outfile = r'C:\folder\arr.asc' #Adjust path
rl = QgsProject.instance().mapLayersByName('DSM')[0] #Adjust to match your raster layer name
e = rl.extent()
h = rl.height()
w = rl.width()
...
0
If you want to use data from a W*S then you should be using a WCS not a WMS. a WMS returns a picture of the data while a WCS returns the actual raster data.
MapServer supports WCS so it should be easy.
1
I use this script:
#!/bin/bash
basename=$(echo "$1" | cut -f 1 -d '.')
mask=${basename}_mask.tif
output=${basename}_edt.tif
nodata=$(gdalinfo $1 | grep "NoData" | cut -d "=" -f 2)
gdal_calc.py --NoDataValue=$2 --calc="A!=${nodata}" --outfile="$mask" -A $1
gdal_calc.py --NoDataValue=$2 --calc="A*B" -...
0
Im suffering from the same problem. Tried to read in raster files with CRS set to EPSG:3067. Same message of ETRS89 datum being discarded arises.
One can add manually the information to the crs-string. I just do not believe that this actually fixes the problem.
#EXAMPLE OF ADDING DATUM INFO TO CRS STRING
library(raster)
file = raster("yourFilePath)
...
1
You are doing this:
rasters <- lapply(paste0(mypath, files), raster)
which makes rasters a list of raster objects. Then you do:
for (i in 1:length(rasters)) {
r <- raster(rasters[[i]])
which is calling raster(...) on a raster object (rasters[[i]]).
That creates a new uninitialised raster.
You probably just want to do r <- rasters[[i]] there.
...
2
An old question, but as I couldn't get the plugin mentioned above to work, here another work-around, tested on QGIS 3.12.3.
Add a legend with limited number of discrete classes (enough to
capture the gradient).
In the composer, add the raster layer and legend.
Under symbol, reduce the height of the symbols (1), deselect the ‘draw stroke for raster symbol’ (...
1
Run a little test. Make a 10x10 random raster:
> r1 = raster(matrix(rnorm(100),10,10))
Feed to ks.test:
> ks.test(r1,"pnorm",alternative="two.sided",exact=NULL)
One-sample Kolmogorov-Smirnov test
data: r1
D = 0.052428, p-value = 0.9463
alternative hypothesis: two-sided
Now normally you'd feed ks.test a numeric vector of ...
0
Try drawing a polygon the same size and shape as the raster (cosmetic layer should be fine) and then interrogate the polygon. For a more accurate answer you could draw a larger polygon and clip it using the raster.
0
At Autodesks Autolisp forum I have added two autolisp functions to import points from ESRI ASCII raster format (.asc file) and to filter out unnecessary points.
Since I received negative points let me explain a bit how and when you can use this option.
Ok this is definitively not the best solution since one can in a first place use various converters to ...
0
Generally, the best way to remove such artifacts is to use Shrink-Expand technique.
QGIS supports this via SAGA/GRASS (eg. https://grass.osgeo.org/grass78/manuals/addons/r.grow.shrink.html). You need to shink your areas by 1 pixel and then expand it back by 1.
After shrink you get rid of small thin areas, after expand back you get initial slightly smoothed ...
0
Here is a new package fasterize to speed up the conversion between polygons (here using sf object) to raster. From the polygon first you need to create normal raster to define desired resolution, and than fill in values from polygon.
Here is whole example (available here):
Read packages:
library(sf)
library(fasterize)
Create polygons:
p1 <- rbind(c(-180,-...
1
EDIT: Changed the function so it can now generate number of slices other than multiples of 4.
I've come up with a solution using the sf package and rotating points around the raster. Here's what a slice looks like.
It works perfectly, you can test it out yourself:
#' Sector Masking
#' @description Slices Raster object into sections
#'
#' @param raster ...
1
The approach that worked for me was very straightforward once I figured it out. There may be a better or different approach, but I'll share the solution I came up with in case it is helpful.
Approach
(1) Use the rasterize function in package raster to convert the points to pixels that match the resolution of the original raster. If only presence/absence is ...
5
Given a data frame of coordinates:
points <- data.frame(x=runif(50), y=runif(50))
and a raster:
rast <- raster(xmn=0, xmx=1, ymn=0, ymx=1, res=0.05)
you can use cellFromXY to find the grid cell associated with each point, and thus the number of unique grid cells.
length(unique(na.omit(cellFromXY(rast, as.matrix(points)))))
2
You can construct polygons like the ones shown in your image by converting your points into an sf object and then using st_buffer. To compute the fraction of each grid cell that is covered by the polygons, you can use exactextractr::coverage_fraction. It's important to first dissolve the buffered points into a single polygon using st_union, so that areas ...
3
I highly recommend you add the -I command which will allow you to load the rasters into QGIS faster. this is a command I've used very succesfully to load a 90GB DEM into postgres and display it in QGIS. also the -t flag will help
raster2pgsql.exe -s 2263 -d -C -I -M -l 2,4,6,8,10,12,14,16,20,26,32,64,132 B:\dem.tif -F -t 100x100 dem | psql -d thedb -U me -p ...
3
That is how rasters are stored in PostgreSQL/PostGIS.
You will have a line with an id, and a raster field. Usually rid and rast. If you tile the rasters (speeds up queries considerably) with -t 128x128 for example, you would have one line per tile.
You can add them to QGIS, using the DB manager. You right click on the layer in the DB manager and choose &...
1
From what I can tell, it is because you have two disconnected polygons in a multipolygon for the clip operation. This seems to produce undesired behavior and is not really supported.
You can use the convex hull to get the boundary area of the shapes:
If you use that to clip, you will maintain the resolution of the original raster in the clipped raster.
...
1
It seems that Wagner VII is not converging which is usually an edge condition - I can get Wagner IV (EPSG:54074) to work either with the raster or with vector data.
1
You can simplify the solution (without a function like def getFeatures(gdf))
import rasterio
out_image, out_transform = rasterio.mask.mask(img, poly.geometry, crop=True)
Same with EarthPy
out_image, out_meta = es.crop_image(img, poly.geometry)
0
Solution was to first us pct2rgb.py to convert images to multiband rasters. Then use gdal_translate to add overviews and compression. Images now overlay correctly in WMS and can be rendered as PNG with transparency using QGIS. Hacky Python script below:
#!/usr/bin/env python
import sys
import os
from osgeo import gdal
import datetime
def my_print(text):
...
0
I was unable to figure this out. Since the scenes I am working with do not have any visible cloud cover in them, I don't need to fret about this. But I would like to know how it could be done. However, one easy way to deal with cloud cover could be adding training sites on visible cloud cover and then masking them out based on the cloud cover class developed ...
3
You can open remote dataset with Rasterio. It's really similar to GDAL like illustrated below
import rasterio
mmap_name = '/vsicurl/https://esgf.nccs.nasa.gov/thredds/fileServer/BioClim/WRE/bio11_equiv_2053.tiff'
with rasterio.open(mmap_name) as src:
print(src.width, src.height)
0
Yes, you can. rasterio is a python wrapper for GDAL, so it should be able to read in anything that the GDAL can.
Here is an example of opening a remote file with rioxarray: https://corteva.github.io/rioxarray/stable/examples/COG.html
import rioxarray
mmap_name = 'https://esgf.nccs.nasa.gov/thredds/fileServer/BioClim/WRE/bio11_equiv_2053.tiff'
dataset = ...
0
Here you can find ArcGIS toolbox for calculating SPI.
This toolbox was made with the same algorithm in the spreadsheet.
https://github.com/Hyun-Woo-Jo/CalculateSPI/releases/tag/v1.0
1
It ended up being a bug that was fixed in version 0.0.27
0
You can do that with standard GDAL easily:
import gdal
url="http://storage.googleapis.com/gcp-public-data-landsat"
f1=f"/vsicurl/{url}/LC08/01/204/030/LC08_L1TP_204030_20200607_20200607_01_RT/" + \
"LC08_L1TP_204030_20200607_20200607_01_RT_B4.TIF"
f2=f"/vsicurl/{url}/LC08/01/204/031/...
1
This works for me. If I select GeoTIFF (Floating Point) and download I get a file called MOD_LSTAD_M_2013-01-01_gs_360x180.FLOAT.TIFF. I can load that into R:
> r = raster("./MOD_LSTAD_M_2013-01-01_gs_360x180.FLOAT.TIFF")
> range(r[])
[1] -12 99999
The 99999 value is the ocean or other missing data, so lets set that to NA:
> r[r[]>...
Top 50 recent answers are included
