texts: QGIS attribute form editor crashes on Spatialite source

texts:
I have collected an attribute table (no geometry) together with other stuff in a Spatialite DB for easy portability. The idea was to hand off just one file and DB for field work.

The Problem: The attribute form created uses values from the attribute table inside the Spatialite DB to determine options for the form. This is handy, because the attributes and option for the form can be added to the table easily (new rows). However, QGIS crashes every time I open the info window on a map. If I use a .csv file instead it works properly. But I loose the dynamic, consolidated properties of the workflow...

As there are really no other File DB options, I am left with .csv, which I would rather avoid.

Any suggestions?


texts:
I have a 3.3GB .osm file and I'm trying to convert it to .map with Osmosis with the following

osmosis --read-xml mexico.osm --mapfile-writer file=mexico.map

but I get the java.lang.OutOfMemoryError, I guess because of the file size.

Is there another way to convert it?


texts: Resampling raster with rectangular pixels

texts:
I have a raster file of deforestation but the pixeles are not square. When I project the raster the information of it says that the pixel size is 0.000275322,-0.000271395. I would like to resample the raster such that I have square pixels that have a size of, for example, 0.0025, -0.0025. I want to make this so that I can match the information of the deforestation raster with a slope raster which cells size is 0.0025, -0.0025.

I have not been able to create rasters with non-square pixels (that might also be a solution).


texts: Im trying to understand how QGIS creates bins of data when using graduated symbols. For example if I create a symbol with 2 classes and the values are 1-2 and 2-3. What will happen with a feature with a value of 2? In what bin will it go?


texts: How to get "all users that made edits on a bbox" in overpass?

texts: I need all the users that made edits on a bbox. I know about the node (user:"username")(bbox)out meta;. What I need is something like this node (user:".*")(bbox)out meta;. But the problem is, it is not showing any result. I don't think you can you regex when ':' operator is used. How do I go about doing this?


texts: Implementing weighted multiple linear regression with GEE

texts:
I contacted gee recently. I have a question about weighted multiple linear regression in spatio-temporal method. I want to make the number of time series images(total is M) as the number of parameters for multiple regression, and make every neighbor pixel(total is N, such as 20*20,N=400) as regression input. Then, different pixels are endowed with different weights.
At present, I know ee.Reducer.linearRegression() can implement multiple linear regression. Instruction is ImageCollection(total is N Images, M Bands per Image).reduce(ee.Reducer.linearRegression(M,1)), but problem is that this method supposed every neighbor pixel has same weight.
How can I make weights of N add to the function?


texts: How to render Mapbox Vector Tiles from OSM data?

texts:
I'd like to render my own vector tiles in the Mapbox Vector Tile format (MVT) format using OpenStreetMap data (ultimately for use with Mapbox's hot Mapbox GL JS client renderer).

Is TileStache + Mapnik + mapnik-vector-tile + PostGIS the way to do this? How does one configure Mapnik to render MVT tiles?


texts: Selecting multiple field values using range in QGIS

texts:
I have a shape of contour lines that I need to filter using the Height Above Sea Level like this:

contour lines from 500 to 1200 and from 2000 to 2500

I can do it in several steps generating a shape using the first condition and another shape using the second, and finally combining them using these expressions:

"ALTURA_SOB" >= 500 AND "ALTURA_SOB" <= 1200
"ALTURA_SOB" >= 2000 AND "ALTURA_SOB" <= 2500

Nevertheless, I want to know if there is a way to do it in 1 step using the 'Select by Expression' tool.


texts: Changing legend spacing to be further away from map border in QGIS?

texts:
I am using QGIS 3.2.3

I am creating a map in which I have added a legend to (see figure below).

My issue is that I am trying to create a legend which has spacing between the two columns. I have achieved this by spacing --> Column space within the legend properties within map compose. However the spacing within the columns has resulted in my items being pushed to the border of the map and legend.

This has created a large amount of spacing in the center of the legend (indicated by the pink box). I need to find a way of creating the spacing to be closer to the border and the items to move further left. seen in the next figure (where the spacing is shared evenly at the sides and the items move closer inwards).

I have also tried using the align button to move its position and that doesn't work either.


texts: Joining Tigerline Census Data 2010 and American FactFinder xls

texts:
I have a shapefile from TIGERLine data from the Census website, as well as a spreadsheet from the American FactFinder. I've edited the excel sheet to a 97-2003 workbook, and am now able to bring it into ArcMap.

Both have fields which I am attempting to join, but they are in different Field Types. These are the GEOID2 from FactFinder and GEOID10 from the TigerLines. I' have been trying to create a new field in the TigerLines file, and intend to populate the fields using the Field Calculator.

I am trying to get the FactFinder attributes into the TIGERlines shapefile, but Double does not work, and the ID Codes are 11 numbers which is too long for Short or Long Integer. They are currently as Text strings, but if I join them like such, everything ends up as NULL.


texts: I am trying to create a map overview of a zoomed in area in Zambia. I am doing this step by step from a very nice tutorial, however, in the 'Map Frame' dropdown list there is nothing to select (See pic)


texts: Rasterize doesn't preserves the class attributes original values

texts:
I'm using this GDAL based rasterize function in python in order to create raster from shapefile based on numerical class values.
I am using:
function of GDALMy class attributes are very long numbers :

so I use those clases for the rasterizing as following (the code is from here: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html):
#Rasterize function
def rasterise_me(raster, vector, attribute,
                fname_out="", format="MEM"):
    """Rasterises a vector dataset by attribute to match a given
    raster dataset. This functions allows for the raster and vector
    to have different projections, and will ensure that the output
    is consistent with the input raster.

    By default, it returns a handle to an open GDAL dataset that you
    can e.g. `ReadAsArray`. If you want to generate a  GTiff on disk,
    set format to `GTiff` and `fname_out` to a sensible filename.

    Parameters
    ----------
    raster: str
        The raster filaname used as input. It will not be overwritten.
    vector: str
        The vector filename
    attribute: str
        The attribute that you want to rasterize. Ideally, this is
        numeric.
    fname_out: str, optional
        The output filename.
    format: str, optional
        The output file format, such as GTiff, or whatever else GDAL
        understands
    """
    # Open input raster file. Need to do this to figure out
    # extent, projection & resolution.
    gdal.UseExceptions()
    g = gdal.Open(raster)
    geoT = g.GetGeoTransform()
    nx, ny = g.RasterXSize, g.RasterYSize
    srs = g.GetProjection()
    min_x = min(geoT[0], geoT[0]+nx*geoT[1])
    max_x = max(geoT[0], geoT[0]+nx*geoT[1])
    min_y = min(geoT[3], geoT[3] + geoT[-1]*ny)
    max_y = max(geoT[3], geoT[3] + geoT[-1]*ny)
    # Reproject vector to match raster file
    vector_tmp = gdal.VectorTranslate("", vector, format="Memory",
                                    dstSRS=srs)
    # Do the magic
    ds_dst= gdal.Rasterize(fname_out, vector_tmp, attribute=attribute,
                        outputSRS=srs, xRes=geoT[1], yRes=geoT[-1],
                        outputBounds=[min_x, min_y, max_x, max_y],
                        format=format, outputType=gdal.GDT_Int32)
    return ds_dst

rasterize=rasterise_me(r"29052019.tif", r"shape.shp", "shape ID",fname_out="raster.tif")

that works and generated a raster, thought the raster seemes to have the same classes but I thought it's becuasethe class is very big number:

after that I concat this data about the pixels to big table but then when I check the classes in the big table I can see that they are wrong- they are all the same negative number:

I don't know where do I lose the numbers values or why.
I have tried to check the dtype of this column but it's int64 before and after the rasterize process.
My end goal: to preserve the original class attributes in the final results


texts: ESRI web scene does not load in Story Map Journal

texts:
The web scene I'm trying to embed into my Journal will only load when I'm logged in (be it while I'm in editing mode or if viewing it from another tab). If opening the story in another web-browser or computer, I see the loading bars for a second or so then they disappear and nothing happens. I've read many posts about scenes not loading within a Story Journal due to HTTPS vs HTTP and shared/unshared issues. Reading over those, I've followed the advice but to no avail.

Here are the details of my organization, Journal and web scene in question, maybe someone can spot the issue:

I've checked that the browsers support WebGL
In my account profile, my organization URL includes HTTPS (https://UTEP.maps.arcgis.com)
The Journal and web scene are both shared publicly.
The Journal URL includes HTTPs as per the it's settings page ( https://utep.maps.arcgis.com/apps/MapJournal/index.html?appid=c2b657d045a04df89fd5911e6178b71c )
The Web scene url includes HTTPs ( https://utep.maps.arcgis.com/home/webscene/viewer.html?webscene=cac10a51b6414913a672559032dd1573&ui=min )
In the configuration of the web scene url for the main stage, I leave "unload when reader navigates away" uncheck and by default, "Load page over a secure connection (HTTPs)" is checked.
The page that is supposed to have the scene embedded is the second after the home page (for those attempting to check it out)
Not sure if this is anything, but ever since I started trying to embedded the web scene and playing around with it, the usage of the underlying scene layer has skyrocketed (I'm talking ~8,000 requests in the past five days).

Is this normal or could it be contributing to the problem?

Is there a request every single time I pan, tilt or zoom around?


texts:
I see that OData V3 has support for geospatial types.

Has anyone published a feed that uses the new geospatial data types?


texts: How to convert lines to polygons in QGIS?

texts:
I have a shapefile with linestring  type which are the roads and I want to be able to use this in PowerBI but  it doesn't support linestring types. So I want to covert each of the lines into polygons. Polygons that are just rectangles of each of the lines as I still want to represent roads.

Is there a way I can do this in QGIS?


texts: Animate layer image changing source (StaticImage)

texts:
I need recreate animated radar (in timeline) with openlayers 3(or 4).
 I try changed source on StaticImage, and work ok!.
But when the image comes out a little bit from the screen, the image does not refresh.
My code:

var extent = [-58.51557, -30.244304187, -63.5245788523, -35.24082];
//var extent = undefined;

var source0 = new ol.source.ImageStatic({
    url: '/static/images/setDatos1.png',
    projection: 'EPSG:4326',
    projection: 'EPSG:4326',
    imageExtent: extent
});
var source1 = new ol.source.ImageStatic({
    url: '/static/images/setDatos2.png',
    projection: 'EPSG:4326',
    imageExtent: extent
});

var sources = [source0, source1]

var radar_rosario_layer = new ol.layer.Image({
    opacity: 1,
    visible: true,
    extent: undefined

});
radar_rosario_layer.setZIndex(9000);

var val = 1;
setInterval(function () {
            if (val != 1) {
                val = 1
            } else {
                val = 0
            }
            radar_rosario_layer.setSource(sources[val]);
        }, 1500);


texts: Why doesn't the split tool return any records in the attribute table?

texts:
I've been trying to split a large land use shapefile using a set of study area polygons as the split feature, but when the function finally finishes running (the land use dataset is quite large), all of the shapefiles created by the split are empty. We also tried using another computer, but found the same result, both using ArcGIS 10 with SP2.

Is there any reason why this is the case? Are there any alternate methodologies which would achieve a similar result to split?

Edit: To be more specific: the split field is the FID in text form and we didn't run into any errors relating to the input fields. The study area polygons completely contain the land use.


texts:
In my project I need the paths in a campus but these paths are non-existent so I'm going to walk around and map (I don't know if it's the correct term) the points to build the corresponding paths.

Is there an easier way to do it?


texts: Prepare raster data for wms using geoserver!

texts:
I have some orthophoto of different resolution (urban 8cm pixel, hill 20cm pixel, mountain 35 cm pixel) which cover an area of 28000 km2. Those are in jpg format, already georeferenced.

I need to build a WMS with this amount of data. I need some best practice based on geoserver for building this service. Server requirements, raster format to use etc...

One other thing while working with this is how can i avoid raster overlaps of different resolution for certain area. Because for the same are could be more than one raster with low resolution or higher resolution. I need to get the higher resolution which is for that particular area.


texts: I have prepared a tool using ModelBuilder. The tool mask the raster data using the selected feature layer. It works perfectly fine in ArcGIS Desktop but when I publish it to the server and uses it for the same analysis I get an error saying that the raster file is invalid. I tried even with another dataset but it also resulted in same error


texts: I am trying to create a radiometric map, I have x, y , z coordinates and a value in excel. I want to create a raster map based on this data. How do I go about doing this. I am using ArcGIS 10 with Target/


texts:
I am trying to edit my legend in tilemill to add a "Note". What I did was essentially copy the   legend-source code to do this. However, I would like the line spacing in the "Note" section to be smaller. I have looked at the CartoCSS manual but am sitll not quite getting what I should do.

Here is my legend code. What code should I write to make the line spacing in the "Notes" paragraph smaller?


texts:
Using ArcGIS 10.4.1 I am trying to build a map of the best state camp grounds to go camping in and also to use a telescope.  I have the camp site shapefile from DNR and used lightpollutionmap.info to build a RGB.tif of light pollution for Michigan Lower Peninsula. I then converted the RGB.tif to colormap because I couldn't figure out what to do with the three bands.

I then converted the colormap to a polygon shapefile and I now have a range from 0(least)-7(Most) for the amount of light pollution for 38k ID's... and from here I am lost.

I want to find what camp grounds have the least total value of light pollution around it, say 50km.  So I created a buffer around each camp site and now want it to register the total value of the converted polygon shapefile behind it(presumably the ones with the lowest total value are in the darkest areas for better telescope usage).  I tried dabbling in merge and intersect but I can't seem to find what I am looking for.


texts: Sharing a QGIS project with geotagged images

texts:
I have a QGIS project file with geotagged images. Is there a simple way to share this with a client?

We have come up with two possible solutions so far, which are not great:

Prepare a QGIS project with relative paths. Zip all the files (shapefiles, qgs file etc) and send it to the client. Instruct the client on how to download and install QGIS in order to view map with images.
Create KML/KMZ files so the client can use Google Earth to view the geotagged images.


texts:
I just picked up R this summer.

I created a nested for loop so that nls() will be performed in 24 raster layers and a new file will be created to save the data.

In this example, I have a raster stack of 48 layers worth 2 years of data. My for loop will perform the calculation in each year which consists of 24 layers, write a new ENVI file and proceed to the next year.

I managed to create a files for i = 1 and i = 2. However, R does not stop after creating new ENVI files. It is loading for forever and the temporary files in my C drive is updating from time to time (C:\Users\CandiceYingWoei\AppData\Local\Temp\Rtmp is empty, but the date modified keeps updating).

I tried to use break to stop the for loop. I think it is not necessary because i will never go over nyears since I already set i in (1:nyears).
    if (i > nyears) { break }

How do I stop R after creating the files?

library(sp)
library(raster)
library(ncdf4)
library(rgdal)

setwd("C:/")
new <- stack("82_83_test.envi")
nlayers <- nlayers(new)
nyears <- nlayers(new)/24

for (i in 1:nyears) {
  for (j in (1+(24*(i-1))):(24*i)) {
     fun<- function(x) {
        # a function that perform nls through each pixel in all raster
        # layers in each year and return several parameters (SOS, EOS, etc).
        return (c(SOS,EOS,LOS,SPUDOY,P_TAmplitude))
        }
     equation<-calc(new[[(1+(24*(i-1))):(24*i)]],fun,forceapply=TRUE)
     date<- 1981+i
     name<-paste("Output_",date,".envi",sep="")
     writeRaster(equation,filename=name,format="ENVI",overwrite=T)
     # this is where the break statement is located
     if (i > nyears) { break }
  }
}


texts: Near Analysis takes much longer in ArcPy?

texts:
I have a feature class of about 12 points, a feature class of hundreds of roads, and need to find the nearest road to each point. So the obvious solution (Near tool) takes about 20 seconds, and the results are fine.

When I try to put it in a script (see below), it takes 90 minutes. Also, even though it doesn't crash/error, the results aren't right (the NEAR_FID and NEAR_DIST for each point are both "-1"). I'm running it within ArcMap as a tool, and the run time is all from the Near analysis (it quickly gets through importing the parameters, and then just thinks for a very long time).

I didn't use any of the optional parameters (radius, record angle, record xy, etc.) when running Near either manually or in the script. The data are the same.

Any ideas what is causing the trouble?

import arcpy
from arcpy import env

pointsFC = arcpy.GetParameterAsText(0)
arcpy.AddMessage("Feature Class to be analyzed: " + pointsFC)

spatialRef = arcpy.Describe(pointsFC).spatialReference

outputWorkspace = arcpy.GetParameterAsText(1)
env.workspace = outputWorkspace
arcpy.AddMessage("Output Workspace " + outputWorkspace)

roadsFC = arcpy.GetParameterAsText(2)
arcpy.AddMessage("Roads Layer: " + roadsFC)

# try to find the nearest road
try:
    in_feature = pointsFC
    near_feature = roadsFC
    arcpy.AddMessage("Finding nearest " + roadsFC + " to " + pointsFC)
    arcpy.Near_analysis(in_feature, near_feature)
except:
    arcpy.AddMessage(">> NEAR didn't work")


texts: GDAL reprojection 4326 to 3857 creates gap

texts:
I have an EPSG:4326 tif covering the world that I need to reproject to EPSG:3857 in order to use it with Mapnik to create a map.

When displaying the original tif and leaving Mapnik doing the reprojection, it is working perfectly. But when using GDAL to reproject beforehand, I get a small gap at the right edge of the tif. On the following screenshot of the tif displayed in a mapping app (like Google Maps), you can see it at the middle of the picture:

Here are the commands I used to create the displayed picture:

gdalwarp \
    -co COMPRESS=LZW -co PREDICTOR=2 \
    -multi \
    -t_srs "EPSG:3857" \
    -r cubicspline \
    -tap \
    -tr 10000 10000 \
    world-4326.tif \
    world-3857-10000.tif

gdaldem color-relief \
    -co COMPRESS=LZW -co PREDICTOR=2 \
    world-3857-10000.tif \
    relief-colors.txt \
    relief-world-3857-10000.tif

Is there a way to remove that gap with GDAL?


texts: Calculating total mileage within shape in Leaflet.draw and/or jsts buffer

texts:
I have a working map that allows users to draw polygons and lines and then draw a buffer around them (using this: https://github.com/skeate/Leaflet.buffer).

The map also has some WFS data on the geoserver side. I was wondering if leaflet has the capability of calculating the total number of miles of road data that's within the shape, or how many sq miles of polygon that are within the shape, etc.

This would serve as a rudimentary buffer analysis. The buffer would be created and then a pop-up window would come up saying "X miles of 'road.data', X sq mi of 'polygon.data' in buffer." or whichever pieces of data are elected to be included in the analysis.

I've search for a bit and I can't find anything like this.


texts:
I'm trying to check if a shapefile has the fast feature count.
I have the following C++ code:

#include <iostream>
#include "ogrsf_frmts.h"

int main()
{
    GDALAllRegister();
    GDALDataset *poDS;
    poDS = (GDALDataset*)GDALOpenEx("C:/Users/root/Desktop/2016/GIS Programming/Plugin_new_new/Shapefiles_Karte/Ortsstrasse/temp/buildings - Kopie (2).shp", GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (poDS == NULL)
    {
        printf("Open failed.\n");
        exit(1);
    }
    if (poDS->TestCapability(OLCFastFeatureCount) == true) {
        printf("true\n");
    } else {
        printf("false\n");
    }

    GDALClose(poDS);
}

The code tells me every time "false" even though I am 100% sure the input shapefile has the fast feature count ability.

So I assume that there is a problem in the function call. Isn't this the way you call the TestCapability function?

I am very new to C++ and have only used GDAL with Python so far ..


texts: How to handle real time GPS data with ArcGIS Javascript API?

texts:
I'm developing a web app using Arcgis Server 9.3.1 and Javascript API

I have a database with gps data from several vehicles and I need to display each vehicle real position in a map. I display this points in a Dynamic Layer, but I need to add an Identify action on hover, so I created a graphics layer, and every 10 seconds, I clear the graphics layer and query the database to get the new coordinates values, then, I create and add each point to the graphics layer, so I can hover over the features and ifentify them.

This is already working in IE and Firefox, but in Chrome this causes a memory overflow.

Does anyone have any suggestions to fix this? Thanks.


texts: Are there any books about the psychology of cartography?

texts:
I have the daunting task to do a seminar about the psychology of cartography. I have zero training in this field and would like to ask if anyone has ever tackled this angle?...I would like to generalize this somewhat and start with the fact that humans think in images, and then maybe narrow it down to symbology, meaning and perception of colour and shape, symbolism and also maybe iconography...but I may be going too deep.

Does anyone know of any good site, literature, YouTube links etc. which could inspire me?


texts: MXD Performance - number of layers vs number of datasets

texts:
My mxd currently has 5 buildings layers accessing 5 different buildings datasets (one for each of 5 counties). Most often there would only be a single county's buildings in the current view, if that matters. Assuming there is nothing else causing speed issues, would my mxd likely draw faster (or worse) if I combined the 5 datasets into a single buildings dataset? And if it does, what would be providing the speed increase? Having less layers to display, or having to access fewer datasets?

EDIT:

I appreciate the input on asking whether the data is local or on a network, what data type it is, etc. However, I am more interested in whether there is a significant cost associated with the number of layers in the mxd (such as the map having to iterate through each layer when drawing), whether there is a significant cost associated with the number of separate datasets accessed, and whether one cost is typically more costly than the other. I'm pretty sure that the cost of accessing multiple shapefiles is greater than the cost of accessing the same data as file gdb feature classes, but how do those compare to the cost of rendering layers? I was hoping that someone knowledgeable of the inner workings of ArcMap/ArcEngine would know how these two variables effect map drawing performance.

*All of the data is stored in file geodatabases on the local machine. Currently the building datasets are in 5 separate file gdbs, but I plan to move them all to the same gdb soon.


texts: Error creating or submitting a task in Google Earth Engine

texts:
I am getting the error "Request payload size exceeds the limit: 4194304 bytes" while I am trying to export certain pixel values clipped to the study area.

I am using the following code:

Export.image.toAsset({
image: ex,
description: 'imageToAssetExample',
assetId: 'exampleExport',
scale: 30,
region: geometry,
});

Link to code: https://code.earthengine.google.com/35db26b2991220ae108059ee71cac39f


texts: Nomalization error "Interpolation of X points outside the convex hull"

texts:
Recently I had to update my version of lidR from 2.3 to 3.0 and since then, the script I had for normalizing my files can not finish the process not even if I indicate opt_stop_early()<-FALSE. The code I am using:
# Function that filters all LAS files
ground_filter = function(chunk)
{
  las = readLAS(chunk)
  if (is.empty(las)) return(NULL)

  ground <- filter_ground(las)

  if (is.empty(ground)) {return(NULL)} # if no ground points are found returns NULL

  lasnorm <- normalize_height(las, tin())
  lasnorm <- filter_poi(lasnorm, buffer == 0) # eliminate buffer area before exporting
  return(lasnorm) # if ground points are found returns normalized data
}

opt_chunk_size(ctg_0) <- 0  # needs to be 0 to process file by file
opt_chunk_buffer(ctg_0) <- 100 # needs to be >0 to avoid edge artifacts
opt_filter(ctg_0) <- "-keep_first -keep_class 1 2 3 4 5" # only read poits classified as 1 2 3 4 5 (noise=class 7)
opt_output_files(ctg_0) <- paste(dz_export, "*_dz", sep = "") # export in directory
opt_stop_early(ctg_0) <- FALSE
opt_laz_compression(ctg_0) <- TRUE  # output has .LAZ format

# Creates a list of files with ground points
output <- catalog_apply(ctg_0, ground_filter)
norm <- catalog(unlist(output))

The process is carried out looking normal with some messages that say  Interpolation of X points outside the convex hull defined by ground points (outside the triangulation) failed and returned NAs.
Chunk 202 of 206 (98.5%): state ⚠
Chunk 203 of 206 (99%): state ⚠
Interpolation of 462 points outside the convex hull defined by ground points (outside the triangulation) failed and returned NAs.
Chunk 204 of 206 (99.5%): state ✗
Interpolation of 462 points outside the convex hull defined by ground points (outside the triangulation) failed and returned NAs.
Error: 462 points were not normalizable. Process aborted.
In addition: There were 50 or more warnings (use warnings() to see the first 50)
Execution halted

The normalization is just the beginning of a whole process and after the normalization is finished the process is stopped. How can I avoid these errors?


texts: How to use actual given area (in hectares) of coal mines dataset to style a point layer

texts:
So what I have is a dataset of mines (coal, bauxite, iron ore etc) in Central India and their lat, long. I also have the values for actual area of each of these mines, as per Environmental clearance documents.

I initially tried looking at Google earth satellite imagery to digitize polygons showing the actual extent of each of these mines. However not all of these mines have begun operations. So satellite imagery does not tell me how big these mines will eventually be and I ended up just making a point layer for the mines

I was hoping to prepare a web map (in cartodb) and one in qgis that allows me to style the point layer according to real area of the mine. Hence if the area of the mine is 100 hectares the point should show as a circle of area 100 hectares and accordingly adjust at different zoom levels. Is there an easy way to do this? one which doesn't involve a lot of code.


texts: Build a table with columns of raster values from multiple raster datasets, using Python, GDAL, or PyQGIS?

texts:
I have a series of raster datasets, clipped to the same extent, with matching cell resolutions, each of which represent land surface temperatures (LST) captured by the MODIS sensor.  Each raster stores LST values for sequential time intervals (daily or 8-day intervals) for a year.

I need to build a table with the following specifications:

each row represents a grid cell
one column with a unique grid cell ID
each subsequent column stores LST values from rasters for each time interval (so up to 365 columns max.)

I realize I could do this easily with ArcPy by creating grid cell centroids from one of the rasters, then use the Multi Value to Point tool to build an attribute table for the resulting point dataset with a field for each LST raster.  However, I am building this tool as a PyQGIS plugin, and so need to do this with FOSS libraries and resources (i.e. GDAL/OGR, PyQGIS, etc.)

At this point, I'm a bit unsure as to the best method for doing this.  Here are the alternative methods I've considered:

Convert all of the rasters to xyz files, then use Python's built-in csv library to merge the xyz files together.  I like this option, but I'm still a Python novice, and I haven't been able to figure out how to iteratively merge csv columns together.
Create a point shapefile from grid cell centroids, then try to incorprate Borys Jurgiel's Point sampling tool plugin to extract the values from all of the remaining rasters.  This tool is written as a plugin however, and would take quite a bit of work to fit into my existing code.
Turn the XYZ value into a 1-D Numpy array, and do some kind of Numpy-related magic? (I'm not very experienced with Numpy at this point...)
Use a RDBMS and SQL to do a relational join between csv files.  I'm familiar with PostGIS and SpatiaLite/SQLite, but this will eventually be part of a QGIS plugin, so I'd raster not force users to have to install an RDBMS just to fulfill this one task in my workflow.

Are there any other suggestions as to how I should attack this?  The first option (joining csv files) seems like the best, but maybe I'm over thinking this problem.

UPDATE
Just to follow up on this issue, I decided to go with option 1.  I convert all of the rasters to xyz file, using gdal2xyz.py (included with standard GDAL installations), then I add a header row and unique id column for each xyz file (stored as a CSV).  The final solution to merge these files together into a single CSV table can be found at StackOverflow, here.


texts: Sentinel-1 HH band not loading in Google Earth Engine

texts:
While importing the VV band in Sentinel-1 datasets from GEE, there is no problem when I use the following lines:

var collection = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filterBounds(dam_basin)
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .select('VV');

But replacing the 'VV' with 'HH' leads to the errors in the subsequent lines which worked with 'VV'. The next lines is just thresholding with values less than -14dB.


texts: Determining raster points inside polygon

texts:
I am trying to perform zonal statistics on a input DEM and a polygon shapefile using ArcPy in Python 2.7. I have done this in Python 3, but want to be able to do it using ArcPy tools in Python 2.7. So far, I successfully got a zonal statistics table, but would like to add another feature.

I would like to create a list of all the raster points that are inside each of the polygons. So for example, one of the polygons has 1631 raster points inside its bounds. I would like to append a column on the table (or some other place, doesn't have to be in the zonal statistics table) that is a list of all these points with the value being the elevation from the DEM.

Figure shows what I would like to get, the list of elevation values for all the raster points in each polygon.

Has anyone done this before or know of a way to do it using zonal statistics in ArcPy for Python 2.7 and ArcGIS 10.6.1?


texts: In ArcGIS, exporting a GRIB file to a different format changes the cell values of the original raster

texts:
I downloaded some climate variables from Copernicus as GRIB format. ArcGIS 10.7.1 allows me to directly add them to my table of contents and it shows me, as usual, some maximum and minimum values. The problem is that when I try to export them to a different format, the new raster (I tried to convert them to Tif and Geodatabase file) presents different values compared to the one originally imported. For instance, if the original file presents max and min values of 86.11 and 37.36, respectively, the converted raster (in any format) presents values that range from 97.64 to 10.46! What am I doing wrong?

Thanks for your time!


texts: For ArcGIS tool can Illustration of Item Description be used as Thumbnail in Contents of Catalog window?

texts:
I've made a Python script tool using ArcGIS Desktop 10.0 SP2 and have accessed its Item Description to Update its Illustration by browsing to a *.png file.  The expected Illustration then appears in Item Description.

However, when I use the Catalog window, to access Thumbnail view of the Toolset in which that tool sits, the Illustration does NOT appear as its Thumbnail.  Am I perhaps missing a step?


texts: GeoWebCache: Where to configure the max number of threads when seeding tiles

texts: Where can I configure the maximum number of threads that are used when a seed tiles request is used? The default seems to be 16


texts: Ordnance Survey coordinates mismatch when converting from Lat/Lon

texts:
I have incorporated the batch Lat/Lon to OSGB option in Excel. The source file was on the main Ordnance Survey website:
https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/coordinates
which is called: coordinate calculations spreadsheet (xls).
everything is fine, but unfortunately, there is a quite serious mismatch between the result, we can get in Excel and the result we can get from i.e. gridreferencefinder.com varies around 60-80 meters, as you can see below:

Does anyone know where this problem might come from?


texts: TileWMS layers with openlayers 3 return 400 error with some zoomlevels

texts:
I have 2 TileWMS layers in my system which work well with some zoom levels and do not work with other zoom levels.

each non-working layer return this error:
400: X,Y values for the tile index were calculated to be {546.0380800062516, -354.8241599711327} which had to be rounded to {546, -355} and exceeds the threshold of 10%. Perhaps the client is using the wrong origin ?

the returned x,y values in the error message of each tile are different.

These two layers were working well with all zoom levels of openalerys 2.
Also, I use the same resolutions array which I used with openlayers2.

other layer types are working well with the same projection and the same resolutions array in openalyers3


texts: Using Select Layer By Attribute to return top 3 records

texts:
I have a feature class in a file geodatabase using ArcGIS 10.1.  I want to create a model that selects the three features with the highest values in one field (Average Price). Basically I want my selection to include the three features with the highest average prices.

How do I accomplish this?


texts: Calculating spatial correlation between two datasets using ArcMap and Spatial Analyst

texts:
I am using ArcMap and the Spatial Analyst extension.
I am trying to find a correlation between my raster datasets. For example, for Net Primary Production (NPP) model I used too many input datasets such as solar radiation, soil, NDVI, etc... So I wanted to find which dataset represents the result NPP map most.
I tried to convert them to ASCII format so that I can analyze them in statistical software like SPSS but when I convert it, the output data .txt size became 1.3 GB which is pretty big for a text file. So I couldn't open it with or without any kind of big text file reader software.
Do you have any ideas about it?


texts: Update feature class from values in csv file

texts:
I'm trying to update a feature class with new data from a csv file. My workflow is as follows:

Make a csv reader from csv file
Identify my data columns using their index position
Use an update cursor to update the values in the shapefile with data from the csv

Here is my code so far:

import csv
import arcpy
csvfile = "C:\\Data\\data.csv"
shpfile = "C:\\Data\\TEST.gdb\\azcounties"
file = open(csvfile,"r")
csvReader = csv.reader(file)
header = csvReader.next()
nameIndex = header.index("NAME")
dateIndex = header.index("DATE")
valueIndex = header.index("VALUE")
dateField = "DATE"
nameField = "NAME_2"
valueField = "VALUE"
with arcpy.da.UpdateCursor(shpfile,(nameField,dateField,valueField)) as cursor:
 for row in cursor:
     for rw in csvReader:
         row[0] = rw[nameIndex]
         row[1] = rw[dateIndex]
         row[2] = rw[valueIndex]
         cursor.updateRow(row)

It seems to almost work. I don't get any errors, but only the first row in the shapefile gets updated.


texts: Find nearest neighbours faster in PostGIS

texts:
I currently have a query that calculates the number of useful amenities within a certain radius of each town, using OSM data:

create table places as
select town.name, place, town.way, (
  select count(*)
  from planet_osm_point pub
  where (pub.amenity in ('pub','cafe','restaurant','supermarket','fuel') or pub.shop in ('bakery','convenience')) and pub.way <-> town.way < 5000
) + (
  select count(*)
  from planet_osm_polygon pub
  where (pub.amenity in('pub','cafe','restaurant','supermarket','fuel') or pub.shop in ('bakery','convenience')) and pub.way <-> town.way < 5000
) pubs
from planet_osm_point town where place in ('hamlet','village','town','city');

It works fine, but is very slow, probably because of the high level algorithm: compare every town to every amenity.

Are there any easy ways to optimise this? (Both tables have indexes already.)

"Explain" returns this:

 Seq Scan on planet_osm_point town  (cost=0.00..260497454.31 rows=6490 width=125)
   Filter: (place = ANY ('{hamlet,village,town,city}'::text[]))
   SubPlan 1
     ->  Aggregate  (cost=17614.09..17614.10 rows=1 width=0)
           ->  Seq Scan on planet_osm_point pub  (cost=0.00..17612.32 rows=710 width=0)
                 Filter: (((amenity = ANY ('{pub,cafe,restaurant,supermarket,fuel}'::text[])) OR (shop = ANY ('{bakery,convenience}'::text[]))) AND (((way)::box <-> (town.way)::box) < 5000::double precision))
   SubPlan 2
     ->  Aggregate  (cost=22522.32..22522.33 rows=1 width=0)
           ->  Seq Scan on planet_osm_polygon pub  (cost=0.00..22521.83 rows=197 width=0)
                 Filter: (((amenity = ANY ('{pub,cafe,restaurant,supermarket,fuel}'::text[])) OR (shop = ANY ('{bakery,convenience}'::text[]))) AND (((way)::box <-> (town.way)::box) < 5000::double precision))

Answer

My new query looks like this:

create table places as
select town.name, place, town.way,
  (select count(*)
  from planet_osm_point p
  where p.amenity = 'pub' and st_dwithin(p.way,town.way,5000) as p1) +
  (select count(*)
  from planet_osm_polygon p
  where p.amenity = 'pub' and st_dwithin(p.way,town.way,5000) as p2) pubs
...


texts: Use OGR and GDAL in Python outside OSGEO4W shell

texts:
I have QGIS installed, therefore it installed the OSGeo4W with all functionalities and that works great.

I added all the environment variables as proposed in this post: Is there a way to use gdal functions from OSGeo4W out of their shell?

But what I'm looking for is to actually be able to use my normal Python shell outside OSGeo4W and :

  import ogr

This is since I'm using Apache Spark to make some analysis I don't want to have two different Python environments.


texts: Mismatch between GPS points and data layer at high latitude

texts:
I am currently working on animal movement, based on GPS animal-attached device, at high latitude (north-east corner of Greenland, between 80° and 82° north).

I observe an important mismatch between some of the data (both shapefiles and rasters)I am using and satellite images. I can summarize this by grouping them in two categories :

When displayed again (for example) sentinel satellite imagery, data from MODIS dataset (MOD11, MOD29 and MODIS TERRA LST) do match the coastline from satellite imagery. MOD11 and MODIS TERRA LST CRS is WGS84, EPSG:4326, whereas MOD29 CRS is EPSG:6931 (EASE Grid 2.0 North). The match with coastline is confirmed by consistent animal GPS location on the coast.

When displaying coastline from the GSHHG dataset (https://www.soest.hawaii.edu/pwessel/gshhg/), MASIE Sea Ice Extent (https://nsidc.org/data/g02186) or MASAM Sea Ice Concentration (https://nsidc.org/data/G10005), coastline DO NOT match with satellite imagery. MASIE and MASAM dataset CRS in proj4string format is +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-80 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs, whereas GSHHG CRS is EPSG:4326.

Mismatch is shown in the middle and bottom screenshot shown here (top one shows matching raster).

I used QGIS to display data and illustrate my problem, but I am pretty sure the problem does not come from "on-the-fly" projection since I encounter similar problem when displaying data with R, and data extraction (raster::extract function) does confirm the mismatch. It neither comes from the different resolutions of the datasets I am using.
Note that this mismatch only appears in this northeast corner of Greenland, for what I've observed.
I guess it might come from a difference in the projection method, since I have no other clue of what it could be.

Please, feel free to ask any further details.


texts:
Is there any way to generate a virtual batch (.vrt) in batch? I imagine the coordinates should be in a .txt file and related to the name of the images that will be cropped.

There are more than 1,000 rasters and it would be impracticable to carry out the process one by one.

Here is an example of the algorithm used:

gdal_translate -of VRT -projwin -58.5 -19.0 -57.0 -18.0 \ image.tif image.vrt


texts: Insert variable to query builder of Arcpy

texts:
I have a variable to input to query builder in feature class to feature class arcpy's tool. It doesn't work even I think it is already correct on my script.

The script is like below:

import arcpy

Var = str(input())

# Local variables:
Source = "D:\GIS\Source.shp"
Export = "D:\GIS\SHP"

# Process: Feature Class to Feature Class
arcpy.FeatureClassToFeatureClass_conversion(Source, Export, "01_Export.shp", "\"Field\" = " + Var + ", "", "")

Could anyone solve this problem?


texts:
I downloaded from internet a shp. File containing only one contour line layer. When I was trying to generate slope analysis using “Slope”
Tool, an error message “One or more dropped items were invalid and will not be added to the control” popped out and I could not proceed from this point on. I googled and approached couple of attempts:

Drag dataset into “input raster”, rather the layer itself;
Convert the dataset to interger data type by using “ArcToolbox / Spatial Analyst Tools / Math / Int” tool;
Rasterize the polylines to raster by using “Polyline to Raster” tool;
Export the shp. File to CAD and attempted to use sketchup plugin to do the slope analysis ——GIS showed the exporting is successful while the CAD file is not showing up in the output folder….
…….

None of them worked ….Running out of ideas…..let me know if this explains my question clearly and if you guys know where could the problem be?


texts: Why does ogr2ogr SHP import error out trying to create a database?

texts:
All,

Can anyone tell me what the error is with OGR2OGR when trying to import a .SHP into MSSQLSPATIAL?

ERROR 1: MSSQL Spatial driver doesn't currently support database creation.
Please create database with the Microsoft SQL Server Client Tools.
MSSQLSpatial driver failed to create MSSQL:server=[]\[];database
=[];uid=[];pwd:[]

The import statement I'm using is this:

ogr2ogr -f MSSQLSpatial "MSSQL:server=[]\[];database=[];uid=[];pwd:[]" C:\[path to file]\police_districts_2877.shp -a_srs EPSG:2877 -lco "GEOM_TYPE=geometry" -lco "GEOM_NAME=geom" -nln "CCD_PoliceDistricts" -progress

As you can see, I'm not trying to create a database, and there is one already created. as per the 'creation issue' cited here: http://www.gdal.org/ogr/drv_mssqlspatial.html

I have followed several posts, including this:
Shapefile to MSSQL with ogr2ogr fails to make connection, to no avail...

EDIT: I'm also getting a connection error, but I can't see why...I think my syntax is correct:

ERROR 1: Unable to initialize connection to the server for MSSQL:server=[]\[];database=[];uid=[];pwd:[],[Microsoft][ODBC SQL Server Driver][SQL Server]Login failed for user '[username]'.

Thanks!

-m


texts: ArcReader 10.1 trouble (loading and printing is taking very long ime)

texts:
We've run into printing problems in ArcReader 10.1. Our trouble is with a particular floodplain feature class located in an SDE database. It simply takes forever to load! Then, when printing, it takes at least 10 minutes. Does anybody else have this issue/know a solution?

Thanks


texts: JTS difference op followed by an intersection op not returning expected results

texts:
I have two overlapping polygons and I want to get the line segment(s) where the edges cross.

So I begin by performing a Difference operation. That gives me the area where A is "outside" of B.

Then I do an Intersection op. I expect that to give me the line segment where the two polygons "touch", but as you can see I get a smaller than expected result that doesn't extend along the entire edge where they touch. You can see in the above shot where the green result is, but here it is isolated from the two polygons.

What am I doing wrong? Is there a better way to go about this?


texts: Def Query to not display feature if it isn't attributed

texts:
I have a bunch of features that I don't want to keep checking if something is attributed.  I I want to query it out to not show the geometry if a column has the default string 'NOT ATTRIBUTED' to it.

I can get it to run "Name = 'NOT ATTRIBUTED'" to work, but Name = NOT 'NOT ATTRIBUTED' doesn't work.

It won't run at all using the 'IS'  It has to use the equals sign.  It's on a SQL server.


texts: Calling clip function in pyQGIS from python console?

texts: Is there a way to call the clip function in QQGIS from the python console? It is found under geoprocessing tools in the vector menu.


texts: Stop IGraphicsContainer delete Refresh Automatically

texts:
Below is my code where I am first collecting graphic elements to delete, and then iterating the list to delete the elements. Issue is that the ArcGIS gets refreshed after every delete graphContainer.DeleteElement(item). How can I prevent it for automated refresh.

  IActiveView activeView = pMxDoc.ActiveView;
  IPageLayout pageLayout = activeView as IPageLayout;
  IGraphicsContainer graphContainer = pageLayout as IGraphicsContainer;
  graphContainer.Reset();
  IElement pElement;
  List<IElement> toDel = new List<IElement>();
  pElement = graphContainer.Next();
  while (pElement != null)
  {
    //if (pElement is ITextElement)
    if (pElement is IMapFrame)
    {
    }
    else
    {
      toDel.Add(pElement);
    }
    pElement = graphContainer.Next();
  }
  foreach (IElement item in toDel)
  {
    graphContainer.DeleteElement(item);
  }


texts: Cutting contour lines over raster (building) areas

texts:
I am looking to cut the contours in ArcMap and /or AutoCAD so they don't pass through buildings or roads. Is there a process or code to do this, so it doesn't have to be done manually?

Also, looking to cut the contour lines so that they are able to display elevation labels with no lines running through them.


texts: Programmatically sorting attribute table, adding field and sequentially numbering field based on that sort?

texts:
Basically, I have a bunch of photos from a field data capture software that I am exporting into Arc. I need to create a map for these photo locations labeled in the the order the photos were take.

I will have to sort on the timestamp, add a field and then sequentially number the field based on the sort.

This is a map that will be created after each inspection is done, so I would like to automate the process as much as possible.


texts:
When I try to compile I get

Error   1   Cannot register assembly "C:\Users\Patrick\Desktop\DynamicBiking\CSharp\bin\Debug\DynamicBiking.dll" - access denied. Please make sure you're running the application as administrator. Access to the registry key 'HKEY_CLASSES_ROOT\DynamicBiking.BikingTrackModeCmd' is denied.  DynamicBiking2010

Trying to compile this one:
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/SampleZips/101efef5-9071-49ac-9ce9-c5a5fe5cc896CSharp.zip

I have administrator priveleges


texts: Getting mapserver url from ArcGIS Server?

texts:
I have been creating some web-map applications using GoogleMaps Javascript API v3.

I've been bringing in tiles from ArcGIS servers which is really easy using the API.

The code below loads a sample dataset from ArcGIS.com.

var url = 'http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer';
var agsType = new  gmaps.ags.MapType(url,{name:'ArcGIS'});
map.mapTypes.set('arcgis', agsType);
map.setMapTypeId('arcgis');

I have some clients that already have a public ArcGIS server online and in place.  When I ask them for the MapServer URL, they don't know where to find it or how to generate it for their datasets.

I don't know anything about ArcGIS server and have not been able to find the appropriate info in ESRI help docs.

Could someone point me to some info about generating an appropriate MapServer URL that I could then pass onto my clients?


texts: How to cache and server pre-rendered vector tiles in pbf format?

texts:
I want to serve vector tiles server in pbf format. For that I used Mapnik-vector-tiles which is bundled by avecado. Using avecado_server i can render pbf tiles on the fly using

  ./avecado_server -b 2 test/xml/vector_imposm_nocollect.xml 7789

But the full planet avecado_server is not sufficient to server on the fly because it takes more time to render each tile.
So I need pre-render tiles using ./avecado and cache it in disk and server those cached tiles with separate server.

Avecado generate tiles in pbf format.

For this I prefer to go with Tile Stache. But I didn't know how to use Tile Stache for this scenario.

Kindly suggest me a best way to cache and serve the pre-rendered tiles.


texts: Cut off overlapping polygon areas at their cutting edge

texts:
I have a polygon layer with many circle polygons in it, which are overlapping each other.

I now would like to cut off the overlapping areas at their cutting edge:

I use QGIS 2.2 and couldn't find a tool to do it, all vector tools are only working if there are two different layers to cut (there is now way to convert it to different layers because I have over 100 polygons, which are overlapping each other different times.)


texts: Get list of cities within a bounding box

texts: I need to find a web service that allows me to reverse-geocode but with a bounding box. In other words, I want to give the coordinates for a bounding box of latitude and longitude, and be provided with a list of cities.


texts:
It seems that support for streaming input is wishy-washy in ogr2ogr. For example, why can I do this:

curl "http://data.usgin.org/arizona/wfs?service=WFS&version=1.1.0&
request=GetFeature&typeName=azgs:activefaults&maxFeatures=10" |
ogr2ogr -f "KML" /vsistdout/ /vsistdin/

... but I can't do this:

curl "http://data.usgin.org/arizona/wfs?service=WFS&version=1.1.0&
request=GetFeature&typeName=azgs:activefaults&maxFeatures=10&
outputformat=application/json" | ogr2ogr -f "KML" /vsistdout/ /vsistdin/

The first request uses the "WFS" driver for the source layer, and the second request uses the "GeoJSON" driver for the source layer. The first one works fine, but the second gives me:

ERROR 1: GeoJSON parsing error: unexpected end of data (at offset 6000)
ERROR 4: Failed to read GeoJSON data
FAILURE:
Unable to open datasource `/vsistdin/' with the following drivers.


texts:
Can anyone help me convert .xyz file to a contours? I have recieved some .xyz files from the local council and am hoping to turn them into a terrain for 3d modelling. I have tried the ASCII to feature class but the points do not appear on my map (geolocated and all).

I have emailed the council and am waiting to hear back whether or not they have .dem files but would rather work out how to do it myself.

I have attached the .xyz file so someone may be able see if the data is in the correct format for importing.

XYZ DATA File


texts: Object error while running arcpy.CreateFileGDB_management() outside of ArcGIS Pro

texts:
I have a script that starts by creating a FileGDB by calling arcpy.CreateFileGDB_management.

It runs without error on my own laptop, but raises an error when my colleagues try to run it outside of ArcGIS Pro.

This is the code I'm trying to run:

import arcpy
 arcpy.env.workspace = "C:\\LOCAL\\Project"
 home = "C:\\LOCAL\\Project"
 GDB = "test.gdb"
 arcpy.CreateFileGDB_management(home, GDB)

But it gives the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\management.py", line 20234, in CreateFileGDB
    raise e
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\management.py", line 20231, in CreateFileGDB
    retval = convertArcObjectToPythonObject(gp.CreateFileGDB_management(*gp_fixargs((out_folder_path, out_name, out_version), True)))
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 498, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: Object: Error in executing tool

When run from the Python window inside ArcPro there's no issue, it only happens when it's run from outside.

I've tried to run it in 2 different ways:

In PyCharm by setting the interpreter to C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe
Inserting the code one-by-one directly in C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe

I've tested this on 2 different devices, both giving the same error.

The complete script works on my own laptop but I can't seem to get it to work on other devices. It feels like I'm overlooking a setting somewhere.


texts:
My data is a dataframe of ALTITUDE in the first column, and monthly temperatures (let's say for a year) in the other columns. I used the following code to create a linear model for each monthly temperature as the response variable and altitude as the explanatory variable.

y <- names(data)
models <- lapply(y, function(x) {
  lm(substitute(i~ALTITUDE, list(i = as.name(x))), data = data)})

Therefore, I got 12 models to predict temperature based on altitude. Then I read my raster file, which was a DEM, and changed its name to use predict function later.

raster <- raster("Data/raster/DEM.tif")
names(raster) <- c("ALTITUDE")

I took one model from the 12 created to perform predict function, like this:

model1 <-(models[[1]])
rastert <- predict(raster, model1)

So, the values of altitude from the DEM (raster variable), were taken and used on the model1 to create a rastert which is the temperature stimated raster.

Now, I want to perform the predict function like this for the 12 models at once.
I was trying to use the following for loop, but it is not working. I am new using geodata on R. Can someone suggest a link to check out about looping with rasters?

rastert <- NULL

for(i in models){
  rastert[[i]] <- predict(raster,i)
}


texts: GeoServer, OpenLayers so slow for WFS! Why?

texts:
I am displaying 100 000 points in a layer via WFS protocol. I need about 5 minutes. However, 100 000 does not seems so much...

My points are stored in PostGIS, in a column defined this way : AddGeometryColumn('records','position', 4326, 'POINT', 2)

My GeoServer server tell me almost instantly this:

outputFormat = text/xml; subtype=gml/3.1.1
    resultType = results

So I think it it OK for reading in PostGIS, and send the response to my browser.

But my Firefox/Chrome are very long to display it, so much that it tells me that the page is unresponsive many times. So maybe it is my JavaScript or my browser who is long to display. Here is my JavaScript with OpenLayers:

var wfs = new OpenLayers.Layer.Vector("Record2", {
                    strategies: [new OpenLayers.Strategy.Fixed()],
                    protocol: new OpenLayers.Protocol.WFS({
                    version: "1.1.0",
                    url: "http://localhost:8080/geoserver/wfs",
                    featurePrefix: "cite",
                    featureType: "record",
                    featureNS: "http://www.openplans.org/cite",
                    srsName: "EPSG:4326",
                        geometryName: ""
                  })
        });

        map.addLayers([world, wfs]);
        map.zoomToMaxExtent();

Do you have an idea of the reason for this?


texts:
I have a problem, the grass 7 python shell can't find my rasters. When I try to use the python shell to set the region I get "region not found". I can set the region in the grass shell (and other tools I use can read the images from the grass shell), but I'm trying to develop a looping python code for a big part of my master's thesis and stop on this problem. I have the same problem on two computers, I'm using grass 7 (because of i.segment) on windows vista

The code I type into the python shell is:
grass.run_command('g.region', region='xKM_430_420.1').

Am I typing the code in the wrong way? Should I use another command to set the region when using python? Or is it maybe some environmental variables?

Any ideas?


texts: Problem with installation of QGIS 2.8 in Mac OS X 10.7

texts:
I had followed the installation process of QGIS for installing in my Mac OS X 10.7 (http://www.kyngchaos.com/software/qgis).

However, I had a problem in the moment when I open the Qgis, that could see in the next image.


texts: How to add points using coordinates to QGIS without having to create a CSV?

texts:
Is there a way that I can add single points to QGIS (v3.4) by inputting coordinates directly into QGIS, rather than creating a CSV and adding them that way?

I know there used to be a plugin that could do this in previous versions, but it doesn't seem available in 3.4.


texts: Is there service for exporting georeferenced OpenStreetMap (OSM) maps in GeoTiff format?

texts: Is there a service available for exporting georeferenced OpenStreetMap (OSM) maps in the GeoTiff format?


texts: Do attribute indexes improve layer symbology performance?

texts:
ArcGIS Desktop & Server 10.7.1; Oracle 18c:
For feature classes that are symbolized by a values in a single field (multiple symbols):
Would a layer draw faster in the map if its symbology field had an attribute index?
I know attribute indexes improve performance for things like joins and the WHERE clause. But I’m not sure if it would apply to symbology as well (maybe the symbology mechanism ‘searches’ for values — just like the WHERE clause does, in which case, an index might help).
(It’s a difficult thing to test, since ArcGIS is pretty fast no matter what I do. But in my future PROD environment, I’ll have hundreds of thousands of records in feature classes. So I’d like to design them properly now in DEV, if I can.)


texts: Using Geodesic Measurements for Circle Radii?

texts:
I am currently developing an OpenLayers mapping site. Measurements can be made using a line tool, and an area tool. Both of these are currently set to calculate geodesic measurements as outlined in the OpenLayers API.

I use geodesic measurements rather than planar measurements as during user testing people questioned the tool's measurements for distances they already knew (such as driving between towns).

A new feature of the site is for a user to be able to draw a circle on the map of a set radius. OpenLayers only allows for drawing circles using planar distances, so when a user measures the circle with the geodesic measurement tool the values do not match. In the image below the circle planar radius is 10km, but the geodesic line measurement for the diameter is 12km.

Clearly this will leave a user (and myself) wondering which is correct.

Looking at this answer it is seems most desktop GIS systems "ignore" this issue, and return planar measurements and distances. So what is the best practice in terms of user interface and accuracy to deal with planar and geodesic measurements?

Update

I found this Google example which illustrates the issue of radii and the Mercator projection:

http://maps.forum.nu/gm_sensitive_circle2.html

The JavaScript code to draw the circle is as follows:

    var lat1 = (PI/180)* center.lat(); // radians
    var lng1 = (PI/180)* center.lng(); // radians

    for (var a = 0 ; a < 361 ; a++ ) {
        var tc = (PI/180)*a;
        var y = asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc));
        var dlng = atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(y));
        var x = ((lng1-dlng+PI) % (2*PI)) - PI ; // MOD function
        var point = new GLatLng(parseFloat(y*(180/PI)),parseFloat(x*(180/PI)));
        circlePoints.push(point);
        bounds.extend(point);
    }

Does this circle take the curvature of the earth into account?

Final Update

Working code posted at http://geographika.co.uk/creating-a-geodesic-circle-in-openlayers


texts: ArcPy Spatial Join, use join geometry not target

texts:
My goal is to create a layer of lines representing the closest street segment to a set of points.
The code:
arcpy.SpatialJoin_analysis(
    target_features=point_req,
    join_features=selected_streets,
    out_feature_class=pt_st_join_output,
    join_operation='JOIN_ONE_TO_ONE',
    join_type='KEEP_COMMON',
    match_option='CLOSEST',
    distance_field_name='CLOSEST')

Works almost perfectly, doing the join correctly but the geometry included in the table is from the target layer, points, rather than the join layer, lines.
I could use a cursor to run through the points layer and build a layer of lines by finding each row's joined ID in the lines layer. I'm hoping to avoid this out of concerns about runtime since I'd have to loop through the lines layer for each point.
It seems like there should be a more elegant way to do this.


texts: Average width of irregular polygons QGIS 2.12.0 Lyon

texts:
I'm interested in examining the average width of streets, that are irregular polygons, for a project. I'm using QGIS 2.12.0 Lyon, and have lines (power cables)  and Polygons (streets) as layers. For further analysis I need to know how wide the streets, which overlap with cables, are on average.
I already found a try of a solution for this topic at How can I calculate the average width of a polygon?, but unfortunately I don't get rid of error Messages when I try to apply Peter Krauss's solution in QGIS and I get told it's better to create a new question instead of asking there . It's quite plain to me, that there's no perfect solution for finding the average width for irregular Polygons, but approximatly would be absolutly enough for me. Especially because my knowledge in programming is actually a lot to bad for this.

If the power cables layer overlaps the street layer crosswise, I can use the length of the new layer which is created by intersection.

(street = blue, cable = black, intersect= orange)

But if the cable don't overlap crosswise but along the street, I need the average width of the street and not a punctual selection at one point because the width can change. Also there could be curves or crossroads in the street, what makes it even worse.

(street = blue, affected part of the street layer =purple , cable = black, intersect= orange)

Futhermore I'm actually searching for an way in which you don't have to do that much manual work like drawing lines, but rather a solution in which you can recall the average width of every street any times.

Perhaps someone can explain all necessary steps so I can use it?


texts: Bad calculations with Colombian projection (EPSG:3116) in GeoTools?

texts:
Solved, reading mail-archive I can understand that there is a conflict between:

gt-epsg-wkt-9.2.jar
gt-epsg-hsql-9.2.jar

I tested removing each one, with 'hsql' I get the same wrong result, but with 'wkt' the result is apparently correct:

POLYGON ((1008607.2669034603 934072.2524703469, 1119664.4273319552 934150.3439472048, 1119791.4166697538 823552.3609040213, 1008616.3988750693 823493.7092892379, 1008607.2669034603 934072.2524703469))

I slightly understand that there is a problem with the coordinate axis, as sugested a comment below, but I tried with both libraries:

 System.setProperty("org.geotools.referencing.forceXY", "true");

The result is still wrong, I suppose that must exist a solution without exclude a library gt-epsg-hsql-9.2.jar.

I tried to project a square with KML coordinates:

POLYGON ((-74 4.00000000013149, -73 4.00000000013149, -73 3.0000000000987583, -74 3.0000000000987583, -74 4.00000000013149))

Into the EPSG:3116, Magna-Sirgas-Bogota coordinates, so the result is:

POLYGON ((-9130927.896477975 2771311.8407411724, -9105351.47332166 2884393.3410715847, -9072295.496257499 2876530.4903936666, -9099848.751098553 2763923.9776131, -9130927.896477975 2771311.8407411724))

So, the result is between Cuba and Miami:

The code that i use is:

CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3116");
MathTransform   transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
OutGeometry = JTS.transform(myGeometry, transform);

Is possible that the code below is wrong? or the library has a bug?, i've see the reported bug:
Swiss Issue

And i think that is something like that, what do you think about this issue?, I'm using the Geotools Api 9.2.


texts: Calculating Population Within Buffer in QGIS

texts:
I am new to QGIS and GIS systems in general. I have a shapefile of zip codes with associated populations inside them. I also have random lat/long points. I am trying to find the number of people within a certain radius of these points.

How would I do that using QGIS?


texts: How to get snapped points on the line with v.clean?

texts:
I followed the answer posted by @Micha in splitting a lines layer using points, because my points were not exactly on the line, I used v.distance and v.patch described in the answer.

Then I tied to run v.clean to get snapped point vector, but I tried every cleaning tool as the following figure

but all results looked like the original line I used in v.distance, not point vector. How to use v.clean to get the snapped points?

I am using GRASS 7.0.5 under WIN10 64bit.

Update#1

I tried the answer posted by @Micha and successfully get those snapped coords and add to the original points attribute table via v.db.addcolumn and v.distance, but as I tried the following command, error occurred,

v.out.ascii -c -r --overwrite input=tcpoint@nl columns=snap_x,snap_y format=standard separator=comma precision=20 | v.in.ascii --o input=- output=snapped_pts x=5 y=6 cat=3 columns="east double,north double,id integer,snap_x double,snap_y double" separator=comma skip=1

And if I type the command in GRASS GUI, another error showed up:


texts:
I am creating an application (C#) and I want to connect it to GRASS, mainly run grass scripts from my application.

This might seem to be a repeated question however I am thinking of the possibility that I will no longer open GRASS explicitly (double click the application and the GUI will pop-out).  Also I know that GRASS can create PNG images but can grass do this too without being ran?  I really need GRASS to be hidden from the view so that it would be neat.  BTW I know that ofcourse GRASS would still be needed to be installed in the PC and not completely remove it.

Are there any available tutorials other than the one on the wiki page? Or a sample code?


texts:
I would like to add columns xmin, xmax, ymin, ymax and I take coordinates from column 'geometry'. I have exemplary dataframe with 1 records. Where is the problem?
Data:
   FID_1        GODLO                                           geometry
0      0  7.101.01.21  POLYGON ((7340000.000 5425000.000, 7340000.000...

I use, but I am getting the error:
TypeError: 'Polygon' object does not support item assignment

import geopandas as gpd

budynki = gpd.read_file('/home/gis/siatki_skorowidzowe/s7_2000/2000s7_2000.shp', rows = slice(0,1), encoding ='cp1250')

xlist = []
ylist = []
for polygon in budynki['geometry']:
    for x, y in polygon.exterior.coords:
        xlist.append(x)
        ylist.append(y)
        for x in xlist:
            minx = float(min(xlist))
            maxx = float(max(xlist))
        for y in ylist:
            miny = float(min(ylist))
            maxy = float(max(ylist))
    polygon['xmin'] = minx
    polygon['xmax'] = maxx
    polygon['ymin'] = miny
    polygon['ymax'] = maxy


texts:
Given a set of polygons of property boundary's where there is a field "Network" which shows us which network a property is in, I need to create a tight boundary for each network -it can be multipart as there will be gaps.

My workflow which mostly works is to buffer by 35m and allow "Dissolve result" and then convert to single part.

This gives me something like the following

Which is close but on the outer edge the dark boundaries (which is the 35m buffer) shouldn't be there, except for the red areas below where it's just internal overlapping boundaries.

I have tried to clip this by the property boundaries in the network but this then re-adds all the small gaps over the roads etc.

I have tried the new boundary tool and the minimum bounding geometry tool but the result is as follows for both.

I also tried polygon-to line-to point and then Delaunay triangulation and finally dissolve all polygons in the network but this is not maintain the boundary or the required gaps as well as the buffer-dissolve method.

Any suggestions on what else I can try?


texts: How do I extract a raster's extent in python?

texts:
I am looking for an easy way to extract a raster's extent in python.

Currently, I am using gdalinfo and then parsing out the corner coordinates.  The problem with this approach is that gdalinfo is not consistent in how it reports the corner coordinates.  I am working with a variety of file types (hdf, tif, img, asc, etc).  I need the extent to be in decimal degrees (or DMS - that I can convert to decimal degrees).

This is being used with a PostgreSQL database so I have access to any PostGIS function, gdal, or ArcGIS, although I would prefer a function that does not involve ArcGIS unless it is fast (I am trying to index a large number of files, frequently).  It can be run on either Windows or Unix.

Thank you,
Kevin


texts: How to overlay a raster using OpenLayers?

texts:
Through the example provided by the examples section of the OpenLayers website I was able to replicate the overlay method of an image. But when I try to add a different raster, now present in the filesystem, I can not get the same result. The format of the raster is in PNG.

I'm trying to do something like this: openaerialmap.org. Where I can overlay my own aerial images above a map (in example is OSM, but I will use Bing Maps in the beggining).

Other similar questions that did not help me:

Using Custom Aerial Raster Imagery with OpenLayers 3 or 4
Raster overlay in OpenLayers
OpenLayers 4 adding a georeferenced image
How to define layer order in Openlayers?

I'm using OpenLayers v5.3

JS code:

import Map from 'ol/Map.js';
import View from 'ol/View.js';
import {getCenter} from 'ol/extent.js';
import {Image as ImageLayer, Tile as TileLayer} from 'ol/layer.js';
import {transform} from 'ol/proj.js';
import Static from 'ol/source/ImageStatic.js';
import OSM from 'ol/source/OSM.js';
import {register} from 'ol/proj/proj4.js';
import proj4 from 'proj4';

var imageExtent = [591338.2230843633878976, -3835585.6850060112774372, 591735.8523469158681110, -3834804.7757325535640121];

var map = new Map({
  layers: [
    new TileLayer({
      source: new OSM()
    }),
    new ImageLayer({
      source: new Static({
        url: "/img/test.png",
        crossOrigin: '',
        projection: 'EPSG:32655',
        imageExtent: imageExtent,
        imageSize: [7953,15619]
      })
    })
  ],
  target: 'map',
  view: new View({
    center: [-10997148, 4569099],
    zoom: 4
  })
});

Extent from raster was extracted from QGIS

The console shows no error and the image seems to load, but I can't find antwhere in the map.


texts: String variable substitution in selection expression

texts:
I'm trying to put a variable into a string used to do a selection. At the moment it's like this:

expression = "'No_Props' > 5 AND %s < 200 AND 'ST_Ori' >60" % metres_field
arcpy.SelectLayerByAttribute_management(streets_table, "NEW_SELECTION", '"' + expression + '"')

Printing the expression gives this:

'No_Props' > 5 AND 'Shape_Le_2' < 200 AND 'ST_Ori' >60

which I thought would work.

I've tried several things. Originally the expression was in a docstring:

"""'No_Props' > 5 AND %s < 200 AND 'ST_Ori' >60""" % metres_field

but that didn't work

I considered using the doc attribute but it only worked on a function (not a variable)

Does anyone know what I am doing wrong, and how to make it work?

The expression looks like it should do the trick when it's printed out.


texts: Creating Fence Diagrams from Electromagnetic (EM) Conductance Survey

texts:
I am looking for some help with creating fence diagrams in GIS for Electromagnetic Conductance surveys. I have points along a line that represent sample locations with three different depths (readings) at each location. My goal is to try to create something like a fence diagram along the sample routes using the conductivity results at each of the sample points.

I have a couple pieces of literature as reference, but they are geared more towards geologic modeling. The first is "Tools and Techniques for 3D Geologic Mapping in ArcScene: Boreholes, Cross Section, and Block Diagrams" http://pubs.usgs.gov/of/2014/1167/pdf/ofr2014-1167_carrell-tools-and-techniques.pdf. The second is "Digital Mapping and 3D Visualization/Modelling of Subsurface Geology UsingArcGIS 9.2 and Well Log Data" https://www.scribd.com/doc/3982961/8/How-to-Create-a-Fence-Diagram.

The literature is somewhat helpful but not exactly what I was looking for.

I will be trying to use ArcMap/ArcScene 10.1 with access to 3D Analyst and Spatial Analyst.


texts: Copy raw rasters referenced by a Mosaic Dataset in ArcGIS Pro 2.2.1

texts: Is there a tool in ArcGIS Pro 2.2.1 that allows me to select and copy the raw underlying rasters referenced by a Mosaic Dataset to a folder?


texts: What to do about QSqlDatabase: QOCI driver not loaded

texts:
QGIS - 3.0.1

I am trying to connect to the Oracle database. I am not an expert in QGIS, pyqgis at all. I took help of following code to connect to Oracle database. Connecting directly to an Oracle Database using PyQGIS

I used the above code using my hostname, port, database name, username and password.

I get following message when I run the code through OsGeo4W shell.

C:\Program Files\QGIS 3.0\bin>python-qgis.bat "C:\Arundhati\Automation Test Scri
pt\Radio Signal\QGIS_Test_Run3.py"
Oracle driver not available
QSqlDatabase: QOCI driver not loaded
QSqlDatabase: available drivers: QOCISPATIAL QOCISPATIAL8 QSPATIALITE QSQLITE QM
YSQL QMYSQL3 QODBC QODBC3 QPSQL QPSQL7
Could not open connection.
QSqlQuery::exec: database not open

How to install the QOCI driver? I couldn't find any help online.


texts: Geoserver won't start after setting up CORS (no jetty installed?)

texts:
I tried to enable CORS by uncommenting the necessary lines in the web.xml file of geoserver. However, I didn't know that I also had to make sure having the correct jetty servlets jar file in my WEB-INF/lib folder (as described in this post). As the CORS didn't work I thought I had to restart tomcat, but since I did that I got an error that jetty couldn't be found:

Exception starting filter cross-origin java.lang.ClassNotFoundException: org.eclipse.jetty.servlets.CrossOriginFilter

Now I tried to find out which jetty version I use to get the correct servlets jar file, but somehow I cannot find any installation of jetty in my geoserver folder. I installed geoserver with the war file and running it on tomcat. And as I understood jetty comes together with geoserver, so it should exists...

What should be my next step now?

EDIT: geoserver version is 2.13.2


texts: Does MySQL has function which will return the points inside polygon?

texts: I have MySQL database with polygons, and also the database with points (column longitude and column latitude), is it possible to return all points inside one polygon?


texts: Merging tracks from multiple GPX files into a single track?

texts:
I'd lke to merge multiple tracks from multiple GPS files into a single track in a single GPX file.

I tried the Merge option in Gps Track Editor and GPSBabel, to no avail: In both cases, I end up with multiple tracks in a single file.

This is what it looks like when uploading the merged file to eg. GPSVisualizer:

Is there a way?


texts: Is Esri Maps for Office available offline? AKA, if I build a beautiful PPT with an interactive map, can I bring this presentation to a client or conference and present to them without having to connect to the Internet?


texts:
I will be making PDFs with maps on them.
Two pages, one a4, one a3 format.

Both will print the same thing, just one is going to be larger. Question - what scale do I print on the PDF? The number.
If the bounding box for the object is the same, the scales should be the same on both images, right? (here I dum out).

In example, The width of the bounding box is 6000m.
I can get that the width of a3 image printed is going to be 0.4025m (and 0.2794m for a4).

So how can I get the scale? Is it just 1:6000?


texts:
Currently I want to use the R geospatial tools within QGIS for carry out several models related to creation of random points inside polygons based on attribute values with the purpose of establishing what is the best R model that depicts poblacional distribution within a determined area.
I have installed currently QGIS 3.14 and R 4.0.2 and I have installed Processinng R Provider plugin but I haven't achieved that It works.
Can you help me to find some solution or other alternatives for accomplish this analysis type with these two software?


texts: What can I use to have aerial imaging delivered via slippy URL format

texts:
We have been given a license to aerial imaging from a local supplier. They are available in:
* 1km x 1km tiles in MGA coordinates in both TIFF and ECW file format
* Mosaiced blocks in MGA coordinates in ECW file format.

Our goal is to have these available on an iPad using the Motion-X GPS HD Application and its custom map & map download features.

However we would need to take these tiles and deliver them from a webserver or tileserver(?) in the slippy URL format (eg: http://.../Z/X/Y.png)

Looking for advise on a webserver or tileserver or whatever setup and/or maybe even some off the shelf software package/server software that will allow us to do this easily.

We do have some other maps & charts provided to us as well but they are in GeoTiff format.

Wondering if you could point me in the right direction?


texts:
I am using GRASS 6.4.3 (installed with osGeo4w), and I have this question：

As I type g.list, a function window will pop out, and all functions and results are fine. But when I type g.list rast, I got this message：

  '%GRASS_PAGER%'  is not recognized as an internal or external command, operable program or batch file.

The photo below, left one is the correct output in popout window as I type in command g.list, and the right one is the message I got in command line as I type in g.list rast：

But I have already add the following paths in my path environment variables：

1. C:\OSGeo4W64\apps\grass\grass-6.4.3\bin
2. C:\OSGeo4W64\bin

How can I solve this problem? I am using GRASS 6.4.3 under windows 8.1 64-bits, please give me any suggestion!


texts:
I am displaying from Geoserver a WMS layer into leaflet in a PNG Format:

    var app_coup = L.tileLayer.wms("http://localhost:8081/geoserver/cite/wms", {
        layers: 'cite:vnm_polbnda_adm3_2014_pdc',
        format: 'image/png',
        transparent: true,
        version: '1.1.0',
        attribution: "test"
        });

After this, I query in JSON for a POP Up of the information of my labels:

/// Popup limite_inter_district
    var owsrootUrl = 'http://localhost:8081/geoserver/cite/wms';
    var defaultParameters = {
        service : 'WFS', //
        version : '2.0',
        request : 'GetFeature',
        transparent: true,
        typeName : 'cite:vnm_polbnda_adm3_2014_pdc',
        outputFormat : 'json',
        format_options : 'callback:getJson',
        SrsName : 'EPSG:4326'
    };
    var parameters = L.Util.extend(defaultParameters);
    var URL = owsrootUrl + L.Util.getParamString(parameters);
    var ajax = $.ajax({
        url : URL,
        dataType : 'json',
        jsonpCallback : 'getJson',
        success : function (response) {
        L.geoJson(response, {
                onEachFeature: function (feature, url) {
                    popupOptions = {maxWidth: 250};
                    url.bindPopup("<b>adm3_name:</b> " + feature.properties.adm3_name
                        ,popupOptions);
                }
            }).addTo(map);
        }
    });

I have the following question:

Why the polygons are displaying in blue when I been setting transparent to 'True'?


texts:
I'm having an issue with my popup in my map made using ArcGIS online.

One of my group layers contain 4 layers. When I turn all of the layers in this group layer off, the second layer's popup still displays when I click on the map. Does anyone know why this is happening?

Also, for some reason, when I have all the layers in this group layer turned on, the second layer's pop up (polygon feature) shows before the 1st layer's popup (a point layer). How do I get the first layer's pop-up to show up first before the second layers?


texts:
I am making an atlas with a fixed scale, and want to include circles with a fixed radius (in scaled kilometers) about the center point in the atlas (see image below). I was able to do this by calculating the appropriate scale using the scale bar and a physical ruler. However, it seems like there should be a more precise way to solve this problem.

Within the expression builder of the circle height and width, the scale of the map can be retrieved with:

map_get(item_variable('map'), 'map_scale')

This returns the scale (30,000), however when I try to scale the target distance (2000 m) by this scale it doesn't match length of the scale bar (also 2 km) -- instead it is about 2/3s of the target radius.


texts: WFS to JSON vector layer error when specifying srsname in query string, no features when specifying value of DefaultCRS in OpenLayers 4?

texts:
I'm following the example from the official OL site for WFS, but I cannot get it to work with a different Feature Server when I add srsname= to the query string as specified in the example:

GET http://fbinter.stadt-berlin.de/fb/wfs/geometry/senstadt/re_spielplatzbestand/?service=WFS&version=2.0.0&request=GetFeature&typenames=fis:re_spielplatzbestand&outputFormat=application/json&srsname=EPSG:3857&bbox=1451195.2005082138,6878797.486467866,1529696.0285570899,6909372.297781937,EPSG:3857

Status Code: 400 Bad Request
Response Body:
<?xml version="1.0" encoding="UTF-8"?>
<ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows/1.1" version="2.0.0">
  <ows:Exception exceptionCode="InvalidParameterValue" locator="service">
    <ows:ExceptionText>parameter SRSNAME invalid</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>

Leaving srsname out produces a response with 0 features:

GET http://fbinter.stadt-berlin.de/fb/wfs/geometry/senstadt/re_spielplatzbestand/?service=WFS&version=2.0.0&request=GetFeature&typenames=fis:re_spielplatzbestand&outputFormat=application/json&bbox=1451195.2005082138,6878797.486467866,1529696.0285570899,6909372.297781937
{
  "type":"FeatureCollection",
  "totalFeatures":0,
  "features":[],
  "crs":{"type":"name","properties":{"name":"EPSG:25833"}}
}

I also tried with srsname=EPSG:25833 and appending .EPSG:25833 for the respective bits in url: function(extent) since it is the default projection that appears in GetCapabilities, but got the same result.

https://jsfiddle.net/bapms31j

Without the bbox filter, the URL for GetFeatures returns 1809 features, so it's definitely not a problem from the Feature Service. The same WFS vector layer works flawlessly in QGIS:

Also tried changing version to 1.0.0 and 1.1.0 without success.


texts: How to convert features of layer from POLYGON to MULTIPOLYGON?

texts:
How can I convert features of a layer from POLYGON to MULTIPOLYGON?
I know how to do that using the postgis function ST_Multi but how could I do to same thing for a layer in QGIS?


texts: Calculating changes in slope angle between cells/areas using ArcGIS Desktop?

texts:
I am doing an estimation of erosion with ArcGIS. I use the UniversalSoilLossEquation (USLE) for this.

Now, i want to calculate the L-factor of this equation. I need the FlowLength for this. My problem begins here:
I calculated the Slope of the DEM --> Slope_DEM
Now i want to do the following: I want to calculate the changes in slope angles between cells and create "NoData" if the change is more than 30%.
Example: If one cell has a slope angle of 50° and the next cell has an angle of 25°, I want to create NoData cause the change of slope angle is more than 30% (here it is 50%). After this the flowlength resets

How can I realize this in ArcGIS?

Here is something I found in literature which describes what I want to do very well:

  "This particular problems is addressed by the cutoff slope angle which
  is a user-input value. The cutoff slope angle is defined as the change
  in slope angle from one cell to the next along the flowdirection. This
  value ranges from 0 to 1 for all areas where the slope angle decreases
  from one cell to the next (if the slope angle increases, there will
  definitely be no deposition). Therefore, the user input value will
  range from 0 – 1 and be dependent upon the amount of sediment carried
  by overland flow. For example, an input value of zero will cause the
  slope length to reset every time there is a decrease in slope angle.
  An input value of one will cause the slope length to never reset. In
  an ideal world, this value would be set by an expert familiar with the
  particular area in question. However, this is not always feasible. The
  literature(Griffin,et al., 1988; Wilson, 1986) suggests references
  that a value closer to 0.5 (slopedecreasing by 50% or greater) is
  appropriate.

(source: http://onlinegeographer.com/slope/hickey_slope_length.pdf)


texts:
I have some dwg files that i need to show as background on my map, but the DWG is not need to be vector layers, i need only to show them with out the need of using.
if i am trying to convert it to DXF so the QGIS read it wrong, its convert everything to point and losing all the data.
example for dwg file:

there is any way to show it like it is, without losing all the shaps when trying to convert to vector layer? maybe show it like raster? then what format its should be?


texts: Saving table from PostgreSQL to MapInfo File?

texts:
I have two tables:
1 with geometry in WGS
2 with geometry in Local System

I tried this:

ogr2ogr -f "MapInfo File" output2.tab PG:"host=localhost user=postgres password=123456 dbname=postgis" -sql "SELECT * from filedata WHERE num=1"

and get files.  However, after I open the files in QGIS, I see that file from the table in WGS coordinates is correct but the file from the table in the local system has a shift.

WGS

Local system

The upper part of the polygons are different but its ok.

I add proj line to file GDAL\projlib\epsg but it did not help.
What am I doing wrong?


texts: I am developing an android app for a usb GPS module of UBLOX. Is there any accuracy information in the raw NMEA data which i can use in my GUI? or do I have to use any algorithm to calculate the accuracy?


texts: Snap points to end points of line within 50cm radius in MapInfo Pro

texts:
Is there a way in MapInfo Pro to automatically move a point dataset to the end points of a line that's within 50cm of a point?

This is an existing dataset so it's not simply turning on snapping. Ideally anything with an error above 50cm should be flagged for human modification.

I have tried the MapCAD tool set but can't get it to do this.

Near tool does this in ArcGIS and there's a partial solution in QGIS in Moving points onto lines (~neighborhood) and Snap points to line and calculate chainage using Open source GIS but I need to use MapInfo Pro for this project.


texts: How come WMS can serve KML? Isn't WMS only meant to serve images?

texts:
If you go to GeoServer, layer preview, by clicking "KML" in the "Common formats" column you get a KML file that contains a URL that can be used to get all information for a layer. For example:

https://chemicalse.etme.gr/geoserver/etme/wms?service=WMS&request=GetMap&version=1.1.1&format=application/vnd.google-earth.kml+xml&layers=etme:LAKES&styles=Lake&height=2048&width=2048&transparent=false&srs=EPSG:4326&format_options=AUTOFIT:true;KMATTR:true;KMPLACEMARK:false;KMSCORE:40;MODE:refresh;SUPEROVERLAY:false&BBOX=20.80842303827691,37.61065467534937,24.69512143309823,39.82801902683197

This is a WMS request (it contains service=WMS). I thought WMS requests only download images, but this request downloads the full data for the layer, including the coordinates of the lines (essentially the layer in vector format). I thought you'd need to use WFS for that.

The thing is that I want to permit a user to view the image of a layer (on a Leaflet map), but I don't want users to be able to download the data for the layer. I deactivated WFS for unauthorized users, but that isn't enough if they can get a KML with WMS.

What have I understood wrong? How can I fix this?


texts: How to toggle group while maintaining sub-toggles

texts:
In QGIS:

I work with large projects in which I have many groups. Often I need to only visualise part of a group. However if I want to toggle off a group and later toggle it on again, all layers in that group are automatically visible. Is there a way to maintain the selection of layers I had visible?

E.g.I want to toggle off this whole group (see image), and when I toggle it on again, to only have the subgroups and layers checked that were checked before.


texts: Open a URL in QGIS from form (or using actions) in Linux

texts:
Does anyone know the correct instructions to open a URL as an action (or from a form) using Arch Linux?

In Windows it's as simple as using 'Open' in the action type, and then [%'field containing URL'%] in the action text.

In Linux this just yields a Malformed URL error. The QGIS guidance suggests xdg-open[%'field containing URL'%], but that just yields Unable to run command xdg-open[%'field containing URL'%]

I've tried the same commands under 'Unix' and 'Generic' action types, but to no different effect.


texts: Can I include polygons that lie across boundaries when using the clip tool?

texts:
I am trying to create a map consisting of two vector layers showing a country divided into districts and contaminated areas within the country.

However, since I'm only interested in the contamination in one specific district, I have used the 'filter' as well as 'clip' tool to narrow it down.

Question: When I clip the layers, all of the contaminated areas that lie across the district boundaries get cut by the boundary line. Is it possible to include all polygons that are either in the district or lie across the boundaries?


texts:
The changes I made to the MXD with ArcPy is not reflecting the changes in ArcMap.

arcpy.env.workspace= "C:\Users\UserID\Desktop"
mxdlocation=os.path.join(arcpy.env.workspace,mxdname)
mxdloc = arcpy.mapping.MapDocument(mxdlocation)

fullPathToNewWorkspace =os.path.join(arcpy.env.workspace,gdb)
for df in arcpy.mapping.ListDataFrames(mxdloc):
  for lyr in arcpy.mapping.ListLayers(mxdloc, '', df):
    if lyr.supports("DATASOURCE"):
        dataName = os.path.basename(lyr.dataSource)
        #using the gdb found within the zip folder.
        lyr.replaceDataSource(fullPathToNewWorkspace, "FILEGDB_WORKSPACE", dataName, True)
arcpy.RefreshTOC()

This is the output of the arcpy codes.

The properties of the layer shows that it belongs to the old datasource


texts: Turf.js Check if polygon "Contains" another polygon

texts:
I am using turf.js for geometrical operation on openlayers3 features but i think it lacks support for "contains" operation.

There are these within and inside operations, but they return "false" when i check it with the polygon. Is there any other way to do this "contains" operation with turf.js or do i have to use the JSTS?


texts: Leaflet: added GeoJSON Feature to FeatureGroup doesn't show up

texts:
I have a geojsonLayer and a selectionLayer. I want to add specific features to the selectionLayer after a click event.

I create the selectionLayer at startup:

selectionLayer = new .FeatureGroup().setStyle(selectionStyle).addTo(map);

To the geojson layer I bind:

layer.on({
    click: featureClick
})

Which contains:

function featureClick(event) {
    event.target.addTo(selectionLayer)
}

This workes though, when I check selectionLayer.getLayers().length, the feature is added. But they doesn't show up on the map at all.


texts: Updating an S-57 base file with the associated update files

texts:
How do I use the update files in the S-57 format (.001, .002, ect.) to update the base file (.000)?
Alternatively, can I see if the base file has already been updated?

Background
My superior supplied me with a collection of S-57 files in a big folder structure which I'm supposed to enter into a PostGIS database. The collection is a large nautical chart that is split into many smaller areas. In the root of the folder structure there's a catalogue file (.031) which I don't know how to read.
Each of the smaller areas has it's own folder structure. Inside those folder structure there are folders numbering from 0 and up. Folder '0' contains the base file (.000) and the X folders contains the update files (.00X). I have no way of knowing if the update files have been applied already and I don't know how to apply them myself either.
I can import the base files into PostGIS without problem with ogr2ogr. I thought I could update the created tables with the -update option, but ogr2ogr cannot open the update files. So I would have to make sure that the base file is updated instead. Which I am currently at a loss for how to do.


texts:
I'm trying to install GDAL osgeo lib into Python 3.6 on Debian Stretch

I installed neccessary packages like:

sudo apt-get update && sudo apt-get install -y \
  binutils \
  libproj-dev \
  gdal-bin \
  libgdal-dev \
  python3-gdal \
  python3-pip \
  python-numpy \
  python-dev \

That installation works for Python 2.7 and Python <= 3.5.

But in interpreter I'm receiving following error:

>>> import osgeo
Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'osgeo'

Do you know any solution for fixing that?


texts:
Using MS SQL Server 2012, when we write a query with

select GEOM.STEnvelope() from table

where table holds approx. 4000 native SQL Server line geometries of type geometry, with 2/6/180 points min/avg/max each, the query takes 18 seconds.

This query is generated by third-party software, so we cannot easily replace it . (Please don't judge the third party devs too harshly, they try for a database-agnostic tool which of course has it's drawbacks.) Ultimately the query tries to get the bbox of the whole table, so the third party tool wraps it in something like

select min(p0.STX) as x1,min(p0.STY) as y1,MAX(p1.STX) as x2,MAX(p1.STY) as y2 from (
  select bbox.STPointN(1) as p0,bbox.STPointN(3)as p1 from (
    select  GEOM.STEnvelope() as bbox from (
      select ... from table
    )
  )
)

The replacement in TSQL (minus some fancy column naming and vertex extraction from the bbox is thus

select geometry::EnvelopeAggregate(GEOM) from table

which, as expected, returns almost immediately.

(How) can we reduce the cost of calling .STEnvelope() on our table?


texts: new administrator - Advice to choise a best web sharing GIS DATA with authorization information level

texts:
i'm a new "GIS 's needer" and i have just Pleasantly surprised to find your community ( thanks google!!!! )
I need for my office's project a appication who permit web sharing(only reading)some gps point on the satellite's map of my country(like google maps)
but also this solution will permit to manage users ...i want just say that there 'll be differents users to read gps map and i want lanage the level of authorization to see information of each point
e.g :when administrator"i" update the map with GPS coordonate of "A" point with some information to represent a location where somebody do a business trip
 I would have wanted that my boss have permission to see all information of the business trip ,but that my office colleague hav just authorization to see the name of man who do business trip...
please sorry for my english i live in french's country and i hope your advice to do my best for this project thank you !!!!


texts: GDAL Java bindings GDAL 2.0.0 in Eclipse

texts:
Does anyone know a quick way to use GDAL Java Bindings in Java Eclipse?

It should be GDAL > 2.0.0 in order to support KEA datatypes. GDAL-Java from OSGEO4W seems to be obsolete.

I am working on Windows 10 on a 64Bit processor and Eclipse Mars.