texts:
I would like to transform tifs in a directory with using python in qgis 3 in order to make tifs with a diffrent crs and adding tfw file.

import glob, os
from PyQt5.QtCore import QFileInfo
from qgis.core import *
import processing
input_f = r'C:\Users\48515\Documents\t1'
output_f = r'C:\Users\48515\Documents\t2'
count  = 0
for raster in glob.glob(input_f + "/*.tif"):

    fileInfo = QFileInfo(raster)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(raster,baseName)

    print(processing.algorithmHelp("gdal:translate"))
    trans = os.path.join(output_f,baseName + "_image_"+str(count)+".tif")
    processing.run('gdal:translate',
                   {'INPUT':raster,
                   'TARGET_CRS':'ESPG:2179',
                   'COPY_SUBDATASETS':0,
                   'DATA_TYPE':0,
                   'EXPAND':'1',
                   'TFW': 1,
                   'OUTPUT':trans})

    count += 1
    iface.addRasterLayer(raster,trans)

However, when i run this code, there aren't any changes. What I quess, there are wrong keys in processing function.


texts: Are latitude and longitude observations independent of each other or correlated?

texts:
I want to develop a bivariate gaussian kernel in order to handle the random noise in my GPS tracks where each point has three parameters (lat, long and time stamp).

I was wondering if latitude and longitude errors are correlated or independent of each other?

My common sense says they are independent because no causal relationship exists between them.


texts: Retrieving projection system used to create shapefile when prj file is missing?

texts:
Possible Duplicate:
  Identifying Coordinate System of Shapefile when Unknown?

Are there any tools to help retrieving the projection system used to create a shapefile when there is no .prj file?

Given that there are well-known projection systems and that most of them deal only with a small area of the earth, it may be possible to find the projection system or just to reduce the possible candidates by providing information about the data extents.  Also, overlaying the data on a map with a known projection system can help to guess the projection system by manually trying them to find one which match.

Does a tool exist to do that?


texts: How to set a default attribute value based on the feature being created

texts:
This question expands on, and updates for QGIS 3, the answers found here.
I have a layer with attributes called "Name" and "Label". I would like to set to default value of "Label" to be the same as that of "Name" unless and until I enter a value for "Label". I used the Attributes Form settings to enter a default value for "Label" of attribute($currentfeature, "Name"); however, since default values are applied only to new features, and the "Name" field has no value until that feature is created, there is no existing value to apply to "Label".
I believe the AutoFields plugin would have worked, as suggested in the linked topic, but this plugin isn't available in QGIS 3. I also tried to use an if()
statement, but again, it refers to attributes of the new feature, so it doesn't set the default for that feature.
As a workaround for my specific task, I know that I can use conditional statements in my label settings to use "Name" if, and only if, "Label" is NULL. So the question here is, how can I set a default value for an attribute in a new feature that refers to another attribute in the same new feature?


texts: How can I fix the error "failed to initialize projection" on Mapnik?

texts:
I have the "failed to initialize projection" problem in the current mac version 3.1.0.
The error text:
RuntimeError: failed to initialize projection with: '&srs900913;'

In older issues you can read: mapnik "use rapidxml instead of libxml2. Only libxml2 supports XML entities" (mapnik/node-mapnik#279 (comment))
One solution is merge all the include files in on file, than the brew version work. But this is very bad!!!
One solution is compile the mapnik self on my mac. But if I do this, how can I link on the compilation process lib libxml2 not rapidxml?
Why is on brew no libxml2 version? Or is anywhere such version?


texts: SQL Server Convert GeoJson format into SQL geometry/spatial data

texts: I'm using javascript and Openlayers 3 and to draw shapes or geometry in osm map.   The data is in geojson format and i want to store it in SQL server 2008 R2 using SQL geometry or spatial data types. Is there a way to convert GeoJSON Format into SQL geometry?


texts: Print preview showing distortion in layers from GeoServer via MapFish and OpenLayers?

texts:
I have an application which shows the print preview before printing the map and the problem is that the layer is distorted in the first image but when I zoom into the map the distortion goes away.

I am using MapFish as the printing module and I assume that OpenLayers deals with the print preview part.

I am new to this so I am not sure if this a code based problem or a setting based problem (either in GeoServer or MapFish)


texts: ArcMap Int tool seems to cut max value in half?

texts:
I have ran into an issue with the ArcMap tools in both 3D Analyst and Spatial Analyst.

I convert an elevation grid from meters to centimeters by using the Times tool with a constant of 100. Works exactly as expected and gives correct values.

When I use the Int tool to convert the floating point numbers to integers, the max value gets cut in half. The max goes from ~23k to 12k.

Here is a picture of the elevation grid at the start (top), after Times (bottom) and then Int (middle):

Any ideas?


texts: ST_Intersects does not use spatial index

texts:
I wrote this part of code:

SELECT CB.id_2 AS CountyID, CA.id_2 AS NeighbourCountyID
 FROM  "Counties10B" CB,"Counties10" CA
WHERE ST_Intersects(CA.geog,CB.geog);

But it does not use spatial index. What could be a problem? Indexes are created on both table.


texts: Using arcpy and script toolbox to add the date to a attribute table

texts:
I would like some information on how to add a date to an attribute table, by adding a user specified date in a script tool that will update the polygons attribute table. In short I have two polygons that I wish to calculate the difference in time, in hours between them. The date is determined by specification inside a script box, seen in the image below.

I would like to use code that goes along with my tool and I need to figure out what I'm missing in the following code; (This is the simple code for the tool)

import arcpy

Mother1 = sys.argv[1]
date1= sys.argv[2] ### I believe date style is in format mm/dd/yyyy hh:mm:ss
Mother2 = sys.argv[3]
date2 = sys.argv[4]

try:
    arcpy.AddField_management(Mother1,"Date","DATE","#","#","#","#","NULLABLE","NON_REQUIRED","#")
    arcpy.AddField_management(Mother2,"Date","DATE","#","#","#","#","NULLABLE","NON_REQUIRED","#")

Then, I tried adding the date to the field such as this;

arcpy.CalculateField_management(Mother1, "Date", date1)
arcpy.CalculateField_management(Mother2, "Date", date2)

But that did not work.

I would like to add the date to the attribute table and then also be able to calculate the difference between the two dates, but need to get the date into the attribute table before I can calculate anything?


texts: What raster smoothing/generalization tools are available?

texts:
I have a DEM that I would like to smooth or generalize to remove topographic extremes (chop off peaks and fill-in valleys). Ideally, I would also like to have control over the radius or level of "blurriness". In the end, I will need a set of rasters that range from slightly blurry to really blurry. (Theoretically, the blurriest would be a constant raster of the arithmetic mean of all values).

Are there any tools or methods I can use (based on Esri, GDAL, GRASS)?

Do I need to home bake my own Gaussian blur routine?

Could I use a low-pass filter (e.g. ArcGIS's filter), and if so, would I need to run it a bunch of times to get the effect of a large radius?


texts: How to query which features of all vector layers intersect with a given bounds?

texts:
I am writing a Openlayers map app, and wondering whether there's a better practice to query which features of which vector layers intersect with a given bounds..

It's certainly not a good idea to iterate features...

thanks for any replies!

UPDATE

In this scenario, vector layers are already added to map, so I assume I don't need to request features from remote GeoServer again... then'll just do intersect check on the fly


texts:
When i try to start GeoDataFarm Plugin in QGIS i get this error :

"2018-11-26T08:32:18     WARNING    Traceback (most recent call last):
              File "C:/Users/robert.banjad/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\geodatafarm\GeoDataFarm.py", line 41, in
              import reportlab
              File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python\qgis\utils.py", line 665, in _import
              mod = _builtin_import(name, globals, locals, fromlist, level)
             ModuleNotFoundError: No module named 'reportlab'

         During handling of the above exception, another exception occurred:

         Traceback (most recent call last):
          File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python\qgis\utils.py", line 337, in startPlugin
          plugins[packageName] = package.classFactory(iface)
          File "C:/Users/robert.banjad/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\geodatafarm\__init__.py", line 34, in classFactory
          from .GeoDataFarm import GeoDataFarm
          File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python\qgis\utils.py", line 665, in _import
          mod = _builtin_import(name, globals, locals, fromlist, level)
          File "C:/Users/robert.banjad/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\geodatafarm\GeoDataFarm.py", line 53, in
          import reportlab
          File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python\qgis\utils.py", line 665, in _import
          mod = _builtin_import(name, globals, locals, fromlist, level)
         ModuleNotFoundError: No module named 'reportlab'
         " and it cannot be loaded, anyone know why?


texts: pgrouting tsp (travelling salesman problem) performance

texts:
I currently use a web API to solve TSP tasks. However my provider limits the query to a maximum of 150 stops. I therefore consider to host my own routing service based on OSM data and pgrouting. My top priority is performance and I was wondering if anybody has experience with the performance of the pgr_tsp function?
I don't care about exact stats or system specs. All I want is a ballpark figure in order to see if pgrouting is an option at all.

E.g. Can I calculate TSP for 1000+ stops and is this a matter of seconds or hours?


texts: Symbology templates mess up colors and legend labels when used in Python script

texts:
I'm having issues when trying to symbolize a layer from a .lyr template using a python script. The templates work when imported in ArcMap as seen below:

But when run with the following Python code:

    df = arcpy.mapping.ListDataFrames(mxd, "dfImage" + str(dfNum))[0]
    legend = arcpy.mapping.ListLayoutElements(mxd, "LEGEND_ELEMENT", "dfImage" + str(dfNum) +"_legend")[0]
    fc1 = arcpy.mapping.Layer(fcGrid)
    arcpy.mapping.UpdateLayer(df, fc1, arcpy.mapping.Layer(symbologyTemplatePath + "VI_data.lyr"))
    arcpy.ApplySymbologyFromLayer_management(fc1, arcpy.mapping.Layer(symbologyTemplatePath + "VI_data.lyr"))
    fc1.symbology.valueField = gridCol
    legend.autoAdd = True
    arcpy.mapping.AddLayer(df, fc1, "BOTTOM")
    arcpy.RefreshActiveView()

    for lyr in legend.listLegendItemLayers():
        legend.updateItem(lyr, styleItem)
    legend.elementHeight = 1.0

It changes the colors, and the legend labels aren't truncated to 2 decimal places:

Does anyone know how to fix this?


texts:
Following my Other Question about non-network integrating routing, I have tried the with-points functions using the code proposed in the answer.
The problem is I don't get the shortest path using these functions
For clarification the pink polyline represents the With-points result and the green one represents Dijkstra's result
Edited
The code I used to  calculate the with_points route is
SELECT *
FROM   pgr_withPoints(
         'SELECT id, source, target, st_length(the_geom)  as cost FROM vr',
         'SELECT pnt.id AS pid,
                 edg.edge_id,
                 edg.fraction
          FROM   points_table AS pnt
          CROSS JOIN LATERAL (
            SELECT id AS edge_id,
                   ST_LineLocatePoint(the_geom, pnt.geom) AS fraction
            FROM   vr
            ORDER BY
                   the_geom <-> pnt.geom
            LIMIT  1
          ) AS edg',
        -888888, -607,
        details := TRUE
      );

vr is my edge table, -888888 and -607 are start and end point id


texts: Cannot use in memory raster object with gdal warp (GDAL Python API)

texts:
These are the operations I am doing in my code:

reading an existing raster (tif)
creating a in memory raster object
reading the existing raster object in chunk as array
making some calculations
writing the array to my in memory object

What I want to do now is, with the in memory raster object ready to be used, feed it to gdal warp, which should also output an intermediate raster object that I'll be use afterwards.

However, when I try to GetRasterBand() out of this latter raster object, I receive an error because it is actually equal to None:

AttributeError: 'NoneType' object has no attribute 'GetRasterBand'

I am pretty sure it is about saving/updating/flushing the status of the raster object after the calculations but I could not find out if this is the actual problem.

How can I use in memory (intermediate) files in my script?

from osgeo import gdal, gdalconst
import numpy as np

in_raster = gdal.Open(r'/my/raster.tif', gdalconst.GA_ReadOnly)
in_b = in_r.GetRasterBand(1)

b1 = in_raster.GetRasterBand(1)
xsize = b1.XSize
ysize = b1.YSize
step = 10  # THIS IS THE CHUNK SIZE!
ystep = ysize / step
ystep = int(ystep)
yresidual = ysize - (ystep * step)

driver = gdal.GetDriverByName("MEM")
temp_calc = driver.Create('',
                          xsize,
                          ysize,
                          1,
                          gdal.GDT_Float32)

temp_b = temp_calc.GetRasterBand(1)

for i in range(step):
    if i != step-1:
        arr_part = in_b.ReadAsArray(
            0, ystep * i, xsize, ystep).astype('float32')

        data = np.nan_to_num(arr_part)

        temp_b.WriteArray(data, 0, ystep * i)
    else:
        arr_part = in_b.ReadAsArray(
            0, ystep * i, xsize, ystep + yresidual).astype('float32')

        # replace nan with zeros
        # https://docs.scipy.org/doc/numpy/reference/generated/numpy.nan_to_num.html
        data = np.nan_to_num(arr_part)

        # write array chunk to (in memory) file
        # https://gdal.org/python/osgeo.gdal.Band-class.html#WriteArray
        temp_b.WriteArray(data, 0, ystep * i)

# get raster resolution
in_gt = in_raster.GetGeoTransform()
in_x = in_gt[1]
in_y = -in_gt[5]

kwargs = {'xRes': in_x,
          'yRes': in_y,
          'targetAlignedPixels': True,
          'resampleAlg': 'average',
          'format': 'MEM'
          }

ds = gdal.Warp('',
               temp_calc,
               **kwargs)
b1 = ds.GetRasterBand(1)
# ds is None and the above line throws the exception


texts: How to measure distance along a line from a specific point

texts: I have shapefiles of a road network and five points.  What I would like to do is measure the distance from each point along the road network using QGIS, but am struggling to achieve this.  So far, I have tried both v.net.distance and v.net.iso, to no avail.  Each attempt returns 'Error:  Vector map not found'.  What could I be doing wrong?


texts: What could I do to make QGIS find all the plugins?

texts:
I'm starting to use Ubuntu 12.10 and would like to get my QGIS running on it. Installation was fine but when I tried to add a package called Points to Paths the Plugin Manager was unabled to find it. Indeed there appeared to be only four plugins available on the two repositories that where activated by default (The QGIS Oficial Repository and the orphans repository). So after activating third part repositories there seem to be more plugins but the desired one doesnt appear.

Is there a way to install these plugins manually on the console? Or what could I do to make QGIS find all the plugins? Has anybody expirienced similar problems?


texts: I wrote a Python script that produces several output features layers in separate .gdb's. As I am testing and modifying my code, I generally go into arcCatalog, find the .gdb under which the feature sits, and delete that feature from the .gdb (leaving the .gdb intact) before rerunning my script. This works fine for all of the feature layers except one. I get a schema lock error because the "file is open in a Python script." I literally have to close my Python IDE (Spyder running from Anaconda) and then go back and delete the feature from ArcCatalog. Any reasons why this is happening and suggestions on how I can resolve it?


texts:
Is there a way to use Qfield the same way as ESRIs Collector? All I want to do i add shapefiles and a tile package of imagery

I'm new to QGIS and Qfield.


texts:
Why does running the script not show the file in the output folder correctly? The script was run in PyQGIS 3.10.
Below find the script used

import processing
import os

path ='DESCUELAJORNADA_CIENT_ALLclip_L8'
output = DESCUELAJORNADA_CIENT_ALLOUTPUT

bandList = [band for band in os.listdir(path) if band[-4:]=='.TIF']

input = [os.path.join(path,band)for band in bandList]

def pca_grass(input,output)
    processing.runAndLoadResults(grass7i.pca,
    {'-f': False,
     'n': False,
     'GRASS_REGION_CELLSIZE_PARAMETER':0,
     'GRASS_REGION_PARAMETER':None,
     'input': input,
     'output':output,
     'percent':99,
     'rescale':[0,400]})
    outList=[band for band in os.listdir(output) if band[-4:]=='.TIF']
    new_output = [ os.path.join(output,band) for band in outList ]
    for band in new_output
        iface.addRasterLayer(band,band)
pca_grass(input,output)


texts: Merge intersecting polygons into one which are part of the same feature

texts:
I'm not a GIS professional so I hope I use the correct wording, at least I hope you understand my question.

Let's say there are five polygons in the feature A. In the shapefile there are
more features. In feature A two polygons from the five existing are intersecting. I want to merge the two so overlapping lines will be deleted and the result is one polygon.

A possible solution would work like this:

I select the feature which 'activate' all five polygons
I switch to edit mode
I select the points of the two polygons
I use function MERGE-SELECTED-POLYGONS

So far I found only functions to merge objects (features?), merge shapefiles, and so on.

Is there a function MERGE-SELECTED-POLYGONS?


texts: inclination / cliff faces from raster heighmap

texts:
I've got a raster heightmap source. I'm working with Quantum GIS.

Extracting contour lines from that is easy, thanks to the plugin. However, I would like to go one step further and find out where on the map the inclination is how steep. More specifically I want to find places where the steepness is above some value x.

Mathematically, I figure calculating the distance between the contour lines would give me what I'm looking for, but is there a plugin or something other that already does this?


texts: How to show number of patches in a neighbourhood?

texts:
Given one raster file containing only one zone that indicates patches of vegetation in a city and a shapefile containing neighbourhoods for the city, explain how you could create a map that shows the number of patches in each neighbourhood.

Do you use topology?


texts:
I am trying to get a simple bar chart on my QGIS 2.6 map. I have done that with diagrams.
It seems to me it is not possible to add simple values to each bar (say there are three bars at a point). I have tried the Dialog plugin, but it only creates a picture file in the print composer. Any ideas?


texts: Densification of Point Layer to get smaller Voronoi Areas

texts:
I'd like to densify a point layer to get a better spatial distribution. With "better" I mean homogeneous distributed. If I had to quantified my improvement, I would calculate the voronoi diagram and after the densification the areas had to be under a special threshold.

Is there an algorithms/tool to do my task. Doesn't matter of ArcGIS, QGIS or Python. I thinks thats a common problem of sensor networks.


texts: st_intersection failing for overlapping multipolygons in sf

texts:
I'm working with some user classification data on a citizen science platform. The way our workflow works is to take multiple user classifications of areas, overlay them on one another, and then look for consensus choice - see https://blog.floatingforests.org/2018/01/05/kelpy-consensus/ for context.

We're transitioning our codebase over to work with sf instead of sp in R for multiple reasons (one being that we've re-written the platform, so we have to re-write the pipeline anyway). As part of that rewrite, we're looking to just stay within the sf ecosystem in order to do out consensus operations instead of going back and forth to rasters (which we could using fasterize, but going from raster to polygon is a PITA). Inspired by this post on gis.stackexchange, I've been trying to use st_intersection, but getting odd results. Here's a reprex.

library(dplyr)
library(purrr)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3

#load the data
test_set <- devtools::source_gist("c95701bd444cda3e342838fd9a090fb3",
 filename = "test_set.R")$value

#visualize
plot(test_set)

So far so good, but when I try and get the intersection...

#get the intersection
overlap <- st_intersection(test_set)

I get the following error:

#> Error in CPL_nary_intersection(x): Evaluation error:
  TopologyException: Input geom 0 is invalid: Self-intersection at or
  near point 314406.13208707399 -5762352.8343122564 at
  314406.13208707399 -5762352.8343122564.

This led me to find the "0 buffer" trick here, but...

#What about if we add a buffer?
test_set_buff <- test_set %>%
  mutate(geometry = map(geometry, ~st_buffer(.x, 0)))

overlap_buff <- st_intersection(test_set_buff)

which throws the following error (even if I up the buffer a bit)

Error in CPL_nary_intersection(x) :
  Evaluation error: TopologyException: found non-noded intersection
 between LINESTRING (312816 -5.76168e+06, 312816 -5.76168e+06) and
  LINESTRING (312816 -5.76168e+06, 312816 -5.76168e+06) at
  312815.57772721699 -5761684.4256249201.

I worried this might be due to invalid polygons, so, I did some digging and pulled up the lovely lwgeom

Sure enough, only two of the multipolygons were above. So, I attempted st_make_valid - which successfully made them all valid - but...

#Test set valid with lwgeom
test_set_valid <- lwgeom::st_make_valid(test_set) %>% st_cast("MULTIPOLYGON")

overlap_valid <- st_intersection(test_set_valid)

Which produces the same non-noded error. I fear that it might be due to some users making selections like

which have some overlap in them - but I'm unsure. It could be something else. Any suggestions for fixes here?

Update After reading this thorough thread, tried the following to both snap to grid and then make valid. Still the non-noded intersection issue.

#Snap to grid at 10m, as Landsat is only a 30m resolution anyway
test_set_snap <- test_set %>%
  mutate(geometry = map(geometry, ~lwgeom::st_snap_to_grid(.x, 10)))%>%
  mutate(geometry = map(geometry, lwgeom::st_make_valid))

overlap_snap <- st_intersection(test_set_snap)


texts: QGIS3 Dark scrollbar customization possible?

texts: In QGIS3, is there any way to make one small GUI customization, which is to make the scrollbar dark on a light background instead of light on a dark background in the Fusion theme?  For some reason I always grab the dark part, must be a habit from using other applications (e.g. Android Studio).  I'm not interested in "Night" mode alternatives, which are too extreme.  The "Windows" theme helps, I go for the scrollbar correctly, but it adds some clutter.  Tweaking the Fusion default theme would be ideal if there is a way to do it, maybe a config file somewhere.


texts: Is kml altitude mode "absolute" equivalent to altitude over WGS84 ellipsoid?

texts:
I am working with an aviation software which outputs points with altitudes over WGS84 ellipsoid.

Now I want to represent those points in Google Earth, writing a kml file.
I know I have to convert the altitudes from ft to m.

Do I have to do any other conversion besides that, or is it "absolute" equivalent to altitude over WGS84?

Documentation is not clear enough for me:

  The absolute altitude mode measures altitude relative to sea level,
  regardless of the actual elevation of the terrain beneath the feature.
  In this way, features can be placed underground, and will not be
  visible. Portions of a feature can extend underground, as in the
  example below. Negative values are accepted, to place features below
  sea level.

  This altitude mode is useful in situations where the altitude value is
  known precisely. GPS tracks, for example, can use the absolute
  altitude mode to display paths created while flying or diving.


texts:
I have created a map of features on an archaeological site. The artifact inventory is in FileMaker Pro 12. Rather than creating a CSV file and using that as the attribute table, I would rather, if possible, use a database as an interface to query the artifact inventory to map the occurrence of items in features to start with and then to work with QGIS analytical tools to see if I can isolate family level information.

QGIS archives said that in 2011 it was not possible, then Albert stated he'd found a way to connect FileMaker Pro to QGIS. I would rather use FileMaker as the interface as it cuts out the CSV state that requires a lot of fiddling.

So my question is can QGIS connect with a FileMaker Pro database through ODBC? If so, what are the steps to link them together so that they may be used? I have enabled ODBC in FileMaker and tried to set up the connection in QGIS but what goes into the blanks on the QGIS pulldown page isn't clear.


texts: PyQGIS - QgsColorRamp Shader, Overwrite discrete values?

texts:
I'm shading a raster layer using a QgsColorRampShader, using discrete values. Basically, I want to make the value '255' a black color because it is a null value. Here is the step-by-step through my code.

First, here is a link to the raster file I'm using: https://github.com/epurpur/PyQGIS-Scripts/blob/master/AestheticMin.tif

In this raster layer, the value range of the pixels is 3 to 900.  First I calculate the value range of the raster layer, because I want to check to know what the values are:

layer = iface.activeLayer()
print("Active Layer: ", layer.name())
provider = layer.dataProvider()
extent = layer.extent()
#Using RasterBandStats to find range of values in raster layer
stats = provider.bandStatistics(1, QgsRasterBandStats.All)
min_val = stats.minimumValue            #minimum pixel value in layer
max_val = stats.maximumValue            #maximum pixel value in layer

value_range = range(int(min_val), int(max_val+1))           #Range of values in raster layer. Without +1 doesn't capture highest value

The important part is the value_range variable, which captures the range of pixel values for the raster layer.

Because each raster will have a different value_range (probably), this is dynamic. Here is how I calculate the first half and second half values for my layer now that I know the value_range:

import numpy

first_half_max = round(numpy.percentile(value_range, 50), 2)
first_half_min = round(min_val, 2)
second_half_max = round(numpy.percentile(value_range, 100), 2)
second_half_min = round((first_half_max + .01), 2)

Now, I shade my raster layer:

raster_shader = QgsColorRampShader()
raster_shader.setColorRampType(QgsColorRampShader.Discrete)

colors_list = [
QgsColorRampShader.ColorRampItem(first_half_max, QColor(204, 219, 255), f"{first_half_min} - {first_half_max}"), \
QgsColorRampShader.ColorRampItem(second_half_max, QColor(153, 184, 255), f"{second_half_min} - {second_half_max}") ]

raster_shader.setColorRampItemList(colors_list)
shader = QgsRasterShader()
shader.setRasterShaderFunction(raster_shader)

renderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(), 1, shader)
layer.setRenderer(renderer)
layer.triggerRepaint()

This produces the following result. This is OK, except there are pixels with a value of '255' (no value) that are lost in this map:

The desired output would look like this.  Notice, the black values around the map and then some pixels inside the map with missing values:

@J.Monticolo's answer below works, but only if the values of the raster layer are less than 255.  If there are values greater than 255, the 255 value is lost. Here is how J.Monticolo creates the color's list (all is the same except the color's list, so I am only including colors_list):

colors_list = [
QgsColorRampShader.ColorRampItem(first_half_max, QColor(255, 255, 255), f"{first_half_min} - {first_half_max}"), \
QgsColorRampShader.ColorRampItem(second_half_max, QColor(153, 184, 255), f"{second_half_min} - {second_half_max}"), \
QgsColorRampShader.ColorRampItem(255, QColor(0, 0, 0), 'No Value') ]

He specifically isolates the value '255' in the last line, where he defines it as a black value. But when using this code with my raster layer, I get the following result:

The legend is correct. But notice, because the value 255 falls within the first_half values, it is shaded white instead of black.

Can I isolate the '255' value while using a range of values to shade the others?


texts: Buffers with Varying Distances around a Set of Points

texts: I have set of points in one shapefile and I have a table with a unique identifier that matches each point. I want to make a buffer around each point, according to the radius that I have specified in the table. Is this possible? I joined the table to the points, but from what I know, the buffer tool doesn't have such an option?


texts: Calculate slope of line segments with QGIS

texts:
I want to calculate the slope of a forest road in mountainous areas.
No problem with ArcGIS following the instructions of the following post:

(see comment of user 'FelixIP')
How do I find the slope of road segments with point elevation data of the same linear path

Summary:
1.) split line (forest road) into equal lengths segments
2.) convert segments to 3D shapes (Interpolate shape tool)
3.) calculate slope of each segment with the field calculator ->

  abs( !Shape!.lastPoint.Z- !Shape!.firstPoint.Z)/ !Shape!.length*100

The output are slope values for every line segment of that polyline. In the attribute table, every segment has its individual slope value now.

This works with ArcGIS. Now I would like to do the same with QGIS.
I tried to use the point sampling tool, tried to rasterize the polyline, tried 'QChainage', the Topology checker and the pProf plugin but it didn't bring the expected results.


texts: Software for amateur topographic surveys

texts:
I'm interested in doing rough topographic surveys using my tripod-mounted Brunton transit compass for azimuth readings and laser rangefinder (and maybe an optical level), so obviously going to be much less accurate than using a theodolite and EDM but find (and fast) for the level of accuracy I need for say 25cm contour maps of rural properties.

I'm looking for low-end (but not necessarily free) Windows GUI software that doesn't have all the bells and whistles of most of the $1,000+ packages out there that I can manually transcribe my data points into and do some sort of least squares network adjustment or whatever it's called to solve traverse misclosures etc as well as be able to easily work with and offset for some of the hacky ways of taking observations like fixed height range rods (painted rebar jammed into the ground) and working in some consumer GPS (3 metre error) readings and hopefully getting a reasonable terrain model out of it.

Suggestions? Also any appropriate reading pitched at this level of topographic surveying and cartography would be welcome too! Like this http://www.lulu.com/shop/stephen-eastmead/landscape-surveying-using-handheld-gps-receivers/ebook/product-20624990.html


texts: Getting HRESULT E_FAIL error when trying to open database connection via connection file

texts:
I developed an Addin in AGP 2.5 where I am using connection files to open database connections at runtime. So far that worked perfect. The connection files are stored on a network drive, referenced in the Addin with absolute paths.
DatabaseConnectionFile connectionFileValidering = new DatabaseConnectionFile(new Uri(AppSettings.Default.sde_validation));
Now that I am doing some maintenance on the Addin I am running into the following error when trying to open the enterprise database:
Database test_db = new Database(connectionFileValidation);
I am getting the following error.
'Database test_db = new Database(connectionFileValidation)' threw an exception of type 'ArcGIS.Core.Data.GeodatabaseException'
Data: {System.Collections.ListDictionaryInternal}
HResult: -2146233088
HelpLink: null
InnerException: {"Error HRESULT E_FAIL has been returned from a call to a COM component."}
Message: "A geodatabase exception has occurred."
Source: "ArcGIS.Core"
StackTrace: " at ArcGIS.Core.Data.Database..ctor(DatabaseConnectionFile databaseConnectionFile)"
TargetSite: {Void .ctor(ArcGIS.Core.Data.DatabaseConnectionFile)}
With respect to the error message I found some stuff that is not related to ArcGIS, but did not help either.
Anybody here can point me the direction of a solution?


texts: L.Control.Search for filtered marker with multiple layer

texts:
I have two geojson layer. tram.js and ubahn.js, both with just marker, I use them in a filter and want to search for these filtered marker. But it isn't working and I don't get any errors.

I used this Plugin: https://github.com/stefanocudini/leaflet-search/blob/master/examples/multiple-layers.html

This is it working: https://labs.easyblog.it/maps/leaflet-search/examples/multiple-layers.html

  var ubahn = {"type": "FeatureCollection","features": [{
      "type": "Feature",
      "properties": {
       "geometry": {
        "type": "Point",
        "coordinates": [
          7.000036,
          51.4452814
        ]
      },
        "name": "Planckstraße",
        "iconcategory": "ubahn"
      }
    },

And

var tram = {"type": "FeatureCollection","features": [{
      "type": "Feature",
      "properties": {
        "name": "Dellwig Bahnhof",
        "iconcategory": "tram"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          6.9251056,
          51.4877009
        ]
      }
    },

And this is my code to filter my marker:

    let checkboxstatus

        var geojsonLayer = L.geoJSON(null, {
        filter: (feature) => {
           const isLineChecked = checkboxstatus.lines.includes(feature.properties.line)
        // CODE OF A WORKING FILTER
},

onEachFeature: function forEachFeature (feature, layer) {
                var popupContent = feature.properties.name;
            if(feature.properties && feature.properties.popupContent) {
                popupContent += feature.properties.popupContent;
                } layer.bindPopup(popupContent);
            },

            pointToLayer: function (feature, latlng) {
                return L.marker(latlng, {
                    icon: L.divIcon({
                        className: feature.properties.iconcategory,
                        html: feature.properties.iconcategory[0].toUpperCase(),
                        iconSize: L.point(25, 25),
                    })
                })

                  .bindTooltip(feature.properties.name);
            }
        }).addTo(map);
function updatecheckboxstatus() {
        checkboxstatus = {
        lines: [],
            // MORE CODE
        }

        for (let input of document.querySelectorAll('input')) {
            if (input.checked) {
                switch (input.className) {
                    case 'line': checkboxstatus.lines.push(input.value); break
                    // AND MORE CODE
                }
            }
        }
    }

    for (let input of document.querySelectorAll('input')) {
        //Listen to 'change' event of all inputs

        input.onchange = (e) => {
            geojsonLayer.clearLayers()
            updatecheckboxstatus()
            geojsonLayer.addData([ubahn, tram])
        }
    }

    /****** INIT ******/
    updatecheckboxstatus()
    geojsonLayer.addData([ubahn, tram]);

And now, I thought I could add the search function and use the geojsonLayer

    /*
    var markerSuche = new L.geoJson([
       geojsonLayer, geojsonOpts
    ]);
    markerSuche.addTo(map); */

    var searchControl = new L.Control.Search({
        layer: geojsonLayer,
        propertyName: 'name',
        container: 'findbox',
        textPlaceholder: 'Suche',
        collapsed: false,
        marker: false,
        moveToLocation: function (latlng) {
            console.coordinates(latlng);
            map.setView(latlng, 12); // set the zoom
        }
    });

    map.addControl(searchControl);

Is it possible, what I had in my head? Did I just miss a ; or ) somewhere? Using the markerSuche layer ( in /* */) isn't working as well. Do I need markerSuche? Where is my mistake? I just want to search for the marker in my filter.


texts: How to define areas from certain NDVI values in QGIS

texts:
I downloaded LANDSAT images from 1985-2018 with Google Engine, 2 for each year. Now I want to know how has the vegetation varied over the years with NDVI, using QGIS 3.6.3. Having calculated the index for an image, I want define areas for certain values of NDVI. I'm trying to get, for example, for a value of NDVI from 0.13 to 0.20 to draw the respective area. How can I do that?
EDIT: NDVI and areas to define.


texts:
I'm looking for an oceans or coastlines shapefile that's a more detailed than what Natural Earth has to offer (1:10m).

I'm working with 2000 Census county shapefiles, but it doesn't cut out lakes or oceans. So the county boundaries are a bit blobby. If I clip it to Natural Earth's oceans, a sliver of New York County gets left inside King's County. (I'm a perfectionist using QGIS and Tilemill).


texts: How to remove trailing zeros from QGIS Graduated Style Class labels?

texts:
I'm sure I've seen a solution to this before but now I can't find it.

When using the graduated fill style in QGIS I get lots of trailing zeros in the labels. This makes my legends look messy in the print composer.

Is there a way to automatically remove these?


texts: Geoserver CSV - changing the default comma separator

texts:
I would like to know if it is possible to change the default comma separator CSV exportation in Geoserver using semicolon i.e. as a separator.

Some of the tables that I published have large text fields which are filled by users. So, these fields can contain commas as a part of data. Once I want to export a table to a CSV file by a WFS request, I am having unstructured content in some rows.

I had a look to the documentation but I've only found that and it seems no useful for my purpose.

Does anybody have any advice to find out a solution?


texts: Rasterio Write Masked Array to TIFF - Shape write and read indices are wrong

texts:
I'm clipping a raster to a shapefile, which returned a masked array, and my goal is to write that masked array to a new .tif file with the correct metadata. I can view the masked array and it looks great. Then, when I write it to a new file using rasterio's write function, something goes wrong, because when I load it back in, it's the wrong shape and it won't plot. See below for details and images:

The original raster and polygon shapefile watershed boundary that I want to clip it to.

Clipped output array and image:

out_img, out_transform = mask(rastDat, shapes=coords, crop=True,indexes=1)

out_mask = np.ma.masked_where(out_img<0,out_img)

fig, ax = plt.subplots(figsize=(8,12))
outPlot = ax.imshow(out_mask,cmap='terrain')

Looks great at this point, and I write it to a new file

#copying metadata from original raster
out_meta = rastDat.meta.copy()

#amending original metadata
out_meta.update({'nodata': -99,
                 'height' : out_img.shape[0],
                 'width' : out_img.shape[1],
                 'transform' : out_transform})

#writing and then re-reading the output data to see if it looks good

with rio.open(outLoc,'w',**out_meta) as dst:
    dst.write(out_fill,1)

newClip = rio.open(outLoc)
clipIM = newClip.read(1,masked=True)

however, the new data is not plotting and has a weird shape. Am I doing something wrong on the writing, or reading, or both? How can I write the raster so that the data will come in and plot like the first figure in this post?

To field any questions, I have tried writing as dst.write(out_fill) as well as reading the data in with and without the mask.


texts: ogrinfo says "Unable to open datasource" with KML?

texts:
I have a KMl file that will open fine in QGIS, but when I try to use ogrinfo I get:

$ ogrinfo -so "myfile.kml"
FAILURE:
Unable to open datasource `myfile.kml' with the following drivers.
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> SDTS
  -> TIGER
  -> S57
  -> DGN
  -> VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> GML
  -> GPX
  -> KML
  -> GeoJSON
  -> GMT
  -> GPKG
  -> SQLite
  -> WAsP
  -> PostgreSQL
  -> PCIDSK
  -> OpenFileGDB
  -> XPlane
  -> AVCBin
  -> AVCE00
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDump
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> PDS
  -> WFS
  -> HTF
  -> AeronavFAA
  -> EDIGEO
  -> GFT
  -> SVG
  -> CouchDB
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> XLS
  -> ElasticSearch
  -> PDF
  -> CartoDB
  -> SXF

Any ideas? I'm working on OSX 10.13,5, if I do ogrinfo --version I get GDAL 1.11.5, released 2016/07/01.

UPDATE: Here's my KML file https://pastebin.com/raw/u6hgrecp


texts: System couldn't find given path while running v.split

texts:
I'm trying to run the Grass Tool v.split and get the following error message:

the grayed outline in German translate to "The System couldn't find the given path"

The .exe-files of the grass Tools are all in the following path:

  C:\Program Files\QGIS 3.4\apps\grass\grass7.6.0\bin

Any idea what problem might be? I get an error message for every gras tool I'm trying to use.

Here is also a Picture of the "Processing Options"

It says "activate; protocol commands; for Vectorlayer use v.external instead of v.in.org; protocol consol output; place of GRASS documentation"


texts: Descriptions for KML Paths on Google Earth for iOS

texts:
I am displaying hurricane paths on Google Earth where if you click on a path it shows a description in a window containing info such as wind speeds, category, location, etc. This works perfectly fine on the desktop version of GE.

My problem comes on the mobile side. Google Earth only displays the path and when I click on a path there is no description. I know this is a long shot since GE hasn't updated their iOS app in ages, but does anyone know a workaround for this or any suggested alternative apps that can resolve this?


texts:
Currently trying to display XY data with what I consider unusual coordinates. These coordinates come from a transit scheduling software called HASTUS. Though I'm able to plot the points in relation to one another it's at a scale of 1:85,222,723,064

Here are a couple of the coordinates from the spreadsheet I was given:

6283072, 1841579
6283413, 1841599

These points are for San Diego where coordinates are around

32.729849, -117.147729

Right now my GCS is WGS 1984 which has worked for me in the past.


texts: How to refresh the HTML, including the linked CSS for a QGIS HTML frame?

texts:
I am using a referenced css file in an HTML Frame in a QGIS 3.4.2 layout.

<link rel="stylesheet" type="text/css" href="map.css">

It works great until I have to go and edit the css.  When I make the change, and press Refresh HTML, the HTML frame does not use the updated css.  If I exit QGIS, and open again...the changes are present.

Is there a way to configure the QGIS project to watch this file? Or manually reload?


texts: Write a query to find all geometries within a distance using ST_DWithin

texts:
I have written a query in PostGIS in which I use the ST_DWithin function to get the geometries which are located in a certain distance:

SELECT images.path, focnetworkpart.__gid
FROM images INNER JOIN focnetworkpart
ON ST_DWithin(images.geom, focnetworkpart.geom, 30000)
WHERE path='x.jpg';

As shown in the query above I use two different tables (images: with point geometry and focneworkpart with line geometry).

There is one more table which has point geometry (nodes table) and I want to include it in the above query.
To be more clear, lets say that I have one input in the image table (with point geometry). What I want is to find all the geometries of node table and focnetworkpart table which are in distance 30000 from the point of the image table.

My question is if I can do this using only one query? Or they have to be two different.

EDITED
This is the table called images (with point geometry)

This is the table called nodes (also with point geometry)

And this is the table called network (with polyline geometries)

What I want is to find all the line and nodes geometries which are within a certain distance from a point in images table.
The output result of the query would be a table with four columns (imageId, PathId, nodeGeometry, PolylineGeometry).


texts: Visualising multiple points using QuickWKT plugin in QGIS?

texts:
Using the QuickWKT plugin in QGIS, a command of the form

POINT (30 10)

draws a point at the desired location. However, how can I achieve to draw multiple points in one go? I am looking for something like this

POINT (30 10, 31 10, 33 10)

to draw three points in one go.


texts: How do I delete data older than 7 days in CartoDB?

texts:
I'm trying to delete data that's older than 7 days. I'm not sure how, though (no SQL experience). Everything I've tried to Google doesn't work.

The dates are formatted like "20200426" if that makes any difference.


texts: Extracting cells from raster by attribute in ERDAS Imagine

texts:
I have a raster file that gives me the slope of all the hills and mountains in an area in degrees (0-90) and I need to extract all the cells which are  >= 45 with an operator in the model builder.  I have a pretty good idea of how this would work, but what I don't know is what operator will extract this data into a new raster file or mask.  Can anyone give me the name of an operator that does this?

Basically, I need to know what operator is ERDAS Imagine's equivalent to ArcMap's Extract by Attributes tool.


texts:
I am new to GIS and looking for the best solution for dealing with anything over 2GB as QGIS struggles on my laptop. I understand it's open source and enjoy using the platform.  But need something with a bit more muscle.

My intended use at the moment is research and analysis.

I'm aware of tools like ESRI/ArcGIS but have never used them.

Are they more suited to larger datasets?

Or will running them from my laptop encounter the same issues, as the real problem is my processor (quad core i7)?

As a Geography student, I qualify for the ArcGIS student pricing which seems reasonable, so I'd be happy to go down that route.

There seems to be so many ESRI product options available that it's difficult to know where to start.

Is anyone able to point me in the right direction?

I've been learning how to use SQL, Python and JS and server-side tech. If custom built options are the best way forward I'm happy to explore that. It'll just take me a lot longer but it's all part of the learning process.

Either way I'm happy getting stuck in!  Just looking for the best way to get started.


texts: More customizable topology options for QGIS

texts:
I'm working on a polyline network where there are a lot of open-ended lines that I do want, but also a lot that I don't. QGIS's topology checker flags all cases, making it a major hassle to go through everything.

Are there better QGIS plugins for topology, or a way that I can narrow my topology rules?

I am currently going through each case manually, but I'm only on error 350 out of 4500 lol

Situation in detail:

Approximately 60k line segments
Every line must be connected to at least one other line
Many lines have only one connection on purpose
Many lines have only one connection because of gaps (mistake)
Lots of actual dangles, but also situations flagged as dangles that are correct
All errors are at a threshold of <1m, but some purposeful lines are also <1m

Edit: I've used GRASS v.clean, and it's better than the topology checker, but still have to do extensive manual checking.


texts: Where can I get zip code shapefiles for years 2000 to 2010?

texts:
I've been searching the web for zip code shapefiles (or equivalents) for years 2000 to 2010, and have only managed to find the following:

http://www.icpsr.umich.edu/TIGER/index.html

which carries the following TIGER/Line releases:
1990 | 2000 | 2000 Urban Area | 2002 | 2004 1st Ed. | 2004 2nd Ed. | 2005 1st Ed. | 2006 2nd Ed.

and http://www.census.gov/geo/www/tiger/shp.html

which carries the following TIGER/Line releases:
2006 2nd Ed. | 2007 | 2008 | 2009 | 2010

These releases contain most of the Zip Code Tabulation Areas that I need. However, the 2002 release in the above link does not appear to contain ZCTAs, and I cannot find releases for years 2001 and 2003.
Does anyone have links to TIGER/Line releases for 2001, 2002, and 2003, or links to other sources that contain zip code shapefiles for those years? I've seen some online vendors that have this data available for a fee, but given the project's budget, free sources would be preferable.


texts: Different results from GDAL raster attribute table command on mac and windows

texts:
I am trying to get the raster attribute table of the following file:
Shift_cult_past_current_prediction_20160107_gr1DD.tif available here: https://www.dropbox.com/s/fmomfpro3j6dmz8/Shift_cult_past_current_prediction_20160107_gr1DD.tif?dl=0

Here is the code I am using:

import gdal
dataset = gdal.Open(dataset_uri) # dataset_uri is path to .tif file
band = dataset.GetRasterBand(1)
rat = band.GetDefaultRAT()

I am using gdal version 2.1.0 on python 2.7.11 (Mac and windows). However, while I get the raster attribute table on windows I get nothing on mac (10.9.5). What could I be doing wrong?

-- EDIT

$ conda list | grep 'gdal'
gdal                      2.1.0               np111py27_1    conda-forge


texts: Re projecting Tiff file from the Mollweide Projection to Mercator Projection using Gdal 2.2

texts:
I'm currently looking at re-projecting a tiff file I have of population density. The file is using The Mollweide Projection.

I have come to understand that Gdal would allow me to do this. I have downloaded the software and have been through a few really good basic tutorials Robert Simmon blog.

I am now trying to apply the same steps to the file I have. I've typed the following command:
gdalwarp -t_srs EPSG:3395 -r lanczos -wo SOURCE_EXTRA=1000 -co COMPRESS=LZW F:\ArcGis\Gdal\NE1_50M_SR_W\DensityMaps\Pop1Warpped.tif F:\ArcGis\Gdal\NE1_50M_SR_W\DensityMaps\Pop1Warpped_Mercator.tif

The cmd window just displays a 0, I've let it run for a couple of hours as I thought it might be a problem of file size. I was thinking that it just meant that the software was processing the file. The original file is 340mb.

I've got a feeling that my gdalwarp command might be wrong or the file size is too big. With my little knowledge of python or gdal the solution is eluding me.

Original File gdalinfo:

C:\Users\ttomt>gdalinfo F:\ArcGis\Gdal\NE1_50M_SR_W\DensityMaps\Pop1Warpped.tif
Driver: GTiff/GeoTIFF
Files: F:\ArcGis\Gdal\NE1_50M_SR_W\DensityMaps\Pop1Warpped.tif
Size is 141988, 60942
Coordinate System is:
PROJCS["World_Mollweide",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mollweide"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",0.0],
    UNIT["Meter",1.0]]
Origin = (-17619594.547443531000000,8750529.461868489200000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-17619594.547, 8750529.462) (  4d35'30.71"W, 83d37'50.88"N)
Lower Left  (-17619594.547,-6484970.538) (107d 3'58.73"E, 55d59' 2.97"S)
Upper Right (17877405.453, 8750529.462) ( 15d11'38.88"E, 83d37'50.88"N)
Lower Right (17877405.453,-6484970.538) (103d21'55.35"W, 55d59' 2.97"S)
Center      (  128905.453, 1132779.462) (  1d17'47.23"E,  9d10'35.76"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray


texts: Fetching qgis plugins on 2.12.0 Lyon is taking a very long time

texts: Opening plugins for the first time on QGIS 2.12.0 Lyon there is a window popping up saying "fetching repositories". You can't do anything until it's finished. It then takes it ages to fetch. Any solutions? (I'm using windows 7)


texts:
I want to add more than one layer to snapping in map

var layerInfos = [

  {layer:results[0].layer}

];

snappingManager.setLayerInfos(layerInfos);


texts: How to add GPInMemoryWorkspace to ArcCatalog?

texts:
I was trying to see if there was a way to visualize the data I am creating and placing in_memory and the help page mentions what the reference is named, but it doesn't at all instruct how to actually reference.

Most literature online talks about interacting with it in ArcPy, but I'm trying to specifically target it in ArcMap or ArcCatalog.


texts:
CMake fails to configure. I get following error:

Could NOT find PythonInterp: Found unsuitable version "1.4", but required
is at least "2.7" (found c:/cygwin/bin/python2.7.exe)

which seems irrational. Probably I have found a cause of this behavior. I think that there is a problem with bat file found in the instruction:

@echo off
set VS90COMNTOOLS=%PROGRAMFILES%\Microsoft Visual Studio 9.0\Common7\Tools\
call "%PROGRAMFILES%\Microsoft Visual Studio 9.0\VC\vcvarsall.bat" x86

set INCLUDE=%INCLUDE%;%PROGRAMFILES%\Microsoft SDKs\Windows\v7.1\include
set LIB=%LIB%;%PROGRAMFILES%\Microsoft SDKs\Windows\v7.1\lib

set OSGEO4W_ROOT=C:\OSGeo4W
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"
path %PATH%;%PROGRAMFILES%\CMake\bin;c:\cygwin\bin

@set GRASS_PREFIX=c:/OSGeo4W/apps/grass/grass-6.4.4
@set INCLUDE=%INCLUDE%;%OSGEO4W_ROOT%\include
@set LIB=%LIB%;%OSGEO4W_ROOT%\lib;%OSGEO4W_ROOT%\lib

@cmd

If I run the bat file and run c:/cygwin/bin/python2.7.exe in cmd I get the following error message:

ImportError: No module named site

Running c:/cygwin/bin/python2.7.exe in cmd without running bat file works fine.

I have found answers to similar questions but none solves my problem. Please help me. What can I change in bat file to make it work?


texts:
Need SQL syntax for ...
I have a table with individual polygons like a jigsaw puzzle.
Given a Lon Lat, I need to find the polygon where the lon lat in located.
I am trying:

select * from "parcels_06003" where "geom" @ ('-120.0077874119654', '38.71168052297362')

geom is my polygon column in the db.


texts: Importing USGS ned30m (NAD83) into Project (NAD83/Conus Albers, OTF on): coords still Lat/Long

texts:
I apologize for posting yet-another-gis-novice's confusion about CRSs, but I have researched the issue and am surprised by what I'm seeing.

In simplest form, I

open a new project
set On-the-fly (OTF) CRS transformation ON
choosing NAD83/Conus Albers (a Projected, and not Geographic CRS) as the Project's CRS (same behavior choosing NAD83/UTM12N)
load a 1x1-degree USGS ned30m raster layer which has an assoc. CRS of NAD83 (Geographic CRS)
note that the layer comes in "skewed" (a North arrow "decoration" points NNE, about 1-o'clock)
note that cursor coordinates appear to be in UTM meters
note that the "ruler" distance-measuring tool appears to measure in meters, showing the raster to be approximately 88 x 112 kilometers
set raster layer's Layer CRS to NAD83/Conus Albers (to match project's)
note that the layer "disappears"
note that coordinates of cursor still appear to be UTM meters
choose Zoom to that layer
note that map is now "squared up" with display (N arrow points N/up)
BEFORE moving cursor over map, note that coordinates still appear to be in UTM meters (in "Coordinate:" box at bottom of display)
AFTER moving cursor over map, note that the coordinates revert to geographic, Lat and Long
after which the "ruler" tool measures the rectangle to be about 1 x 1 (must be degrees again)

As far as I can tell, I have done everything I can to alert QGIS 2.8.6 (similar behavior in 2.14 caused me to try the LTR) to the fact that I prefer meters.

I have also (following other advice I found here) tried Layer->SaveAs (to a Projected CRS), but the resulting 23+GB output hung-up QGIS, causing me to have to abort. (A separate problem owing to my naivete regarding choices I made--or rather, defaults I accepted--in the process of trying to save it.)

This is about as simple an example as I can imagine of wanting, say, to import an elevation file from the USGS and measure the distance around a prominent feature, like a mountain peak, in meters.

Again, I regret my need to broach the subject--again--but I would be deeply appreciative to either find out what I did wrong, or why the behavior I see makes Perfect Sense.  Or both.

--
Oertls

Update: The problem is still present in 2.18.23.
Project->Properties->CoordDisp is set to MapUnits (automatic), but coords are still displaying in Lat/Long, even though the project CRS is Conus/Albers and the NED layer has been saved-out/reloaded--after relaunching QGIS 2.18--as a Conus/Albers CRS.

However IF I set the NED layer CRS to "NAD83" (its default), coordinates display in meters.


texts:
I have downloaded a csv file with grid locations with lon/lat for European Extreme Windstorms footprints (attached). However the coordinate system of this dataset is not given.

As per the data description on the website (http://www.europeanwindstorms.org/database/dataDesc/):

The footprints were created using a model with a rotated pole at longitude = 177.5 degrees, latitude = 37.5 degrees (meaning that the z-axis of the new coordinate system lies at longitude = 177.5, latitude = 37.5 degrees. The footprints are on a regular grid in the rotated coordinate system, with horizontal grid spacing 0.22 degrees.

How to build custom coordinate system in ArcGIS using the information above?

I believe my data should fit the square in the image below:


texts:
I'm using shapely to check if a point is in a polygon, but I also need to check if a point is in a circle. Checking the docs, I don't see how do that...

Maybe I could create a circle?

circle = (x,y).buffer(R)

But it doesn't seem to work using R in miles distance to set the ratio...

Is it possible? Or is there another lib to get it?


texts:
I am working on Modis products. They usually have names like:
 "MYD11A2.A2002185.h20v05.005.2007221163754.hdf".
I am converting these hdf files to geotiff. There is no problem with it. But I also want to build  model that shortens that names such as "2002185.tif" by removing some characters from the beginning and from the end. When I use %name% it gives the same name. Is there anybody having a suggestion?


texts: Calculating distance by common ID (one to multiple) in two point layers?

texts:
In point layer A, the locations of farmers' houses are recorded. Each farmer has a unique ID.

In point layer B, the centroids of farmers' plots are included. Each farmer has more than one plot, thus multiple points in the layer B correspond to layer A by farmer ID.

Based on the common farmer ID, I want to calculate the exact distances (not the nearest) from points of each farmer's house in layer A to his or her multiple plots (points) in layer B using QGIS.

I have tried "Distance matrix" as suggested in most of the posts but it doesn't work well.

Any pointers?

The graph above shows a fraction of my dataset. Green points as layer A for farmer's house; Brown points as layer B for plots. One cluster (visually) can be considered as a village. Usually, plots are clustered around the farmers' house (but not always the case).


texts: Make different buffer around buildings per region - ST_Buffer and ST_Intersection slow query

texts:
I am developing Python code to explore suitable sites for wind power development in some countries using spatial data and Postgresql/PostGIS.

For that reason I have tried to process a buffer around buildings using different buffer values ( 400m - 1000 m) stored in the another table called adminlvl2b. I used OSM data for Germany and Portugal to test the performance.
In case of Portugal the query for one region lasts 1.30 minutes, but for 20 regions the process lasts
1 hours and 40 minutes.  In case of Germany I am not able to process the buffer due to higher number of features stored in the memory while processing ST_Union. So I got an error.

I also tried ST_Collect instead of ST_Union, but it took also 1hour 40 minutes and produced only 3 buffered building in 3 regions although it supposed to produce 20 different buffered buildings per region.
I am aware that St_intersection requires much of the total time.

So my question is, how I can rewrite the query to speed up the process and it will work regardless the size of the building data.

I am working with PostgreSQL 9.5.6 on x86_64-pc-linux-gnu, compiled by gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609, 64-bit

Here is the query used and a map representing the results for Portugal:

DROP TABLE IF EXISTS buffered.buildings_residential;

CREATE TABLE buffered.buildings_residential (pkey serial, geom geometry);

INSERT INTO buffered.buildings_residential (geom)

SELECT ST_Union(ST_Intersection(ST_Buffer(m.geom, a.buffer_m), a.geom))
  FROM selected.buildings_residential as m, selected.adminlvl2b as a

WHERE
  ST_Intersects(m.geom, a.geom)

--and   a.pkey=19    -- in case of one region
GROUP BY a.pkey;
ALTER TABLE buffered.buildings_residential ADD CONSTRAINT buildings_residential_pkey PRIMARY KEY (pkey);

CREATE INDEX ON buffered.buildings_residential USING gist (geom);

ANALYZE buffered.buildings_residential;


texts: Carrying over attributes after building geometry in ArcPy with ArcGIS Pro

texts:
I've created polygons from a table that has two columns- an OBJECTID and a string of coordinates (Polygon_Coordinates). The resulting polygons (called "testing" in the code below) are in the correct locations, but I can't seem to figure out how to carry over the ID, which I need for a join later in the workflow.
Any idea on how to build geometry, but maintain the attributes?
fields = ['Polygon_Coordinates', 'OBJECTID']

with arcpy.da.SearchCursor(table, fields) as table_cursor:
    for row in table_cursor:
        #skip null rows
        if row[0] != '-9999':

            #turn string of coordinates into a list but split them at the commas
            coords= list(row[0].split(','))

            #break up list into list of pairs
            pairs = [coords[i:i + 2] for i in range(0, len(coords), 2)]
            coord_ints = [[float(i) for i in p] for p in pairs]

            #create insert cursor for geometry creation
            geocursor = arcpy.da.InsertCursor(testing, ["SHAPE@"])
            array = arcpy.Array()

            #grab pairs from list, put into array
            for a in pairs:
                array.add(arcpy.Point(a[1],a[0], ID=row[1]))

            #set crs, build polygon
            wgs = arcpy.SpatialReference(4326)
            poly = arcpy.Polygon(array, wgs)
            geocursor.insertRow([poly])


texts: Mapserver - Dynamically Change OGC WFS Filter from WMS BBOX?

texts:
Is it possible to dynamically change the content of an OGC WFS filter based on the BBOX of an incoming WMS request?  I am attempting to serve the response of a WFS which require that I use an OGC filter on each request for the desired result.

A consequence of this is that the BBOX of the WFS request must be contained within the filter as well.

Filtering client side is not possible in this case because 1) the dataset is very large and 2) the filter parameters materially affect the ouptut.

I apologize that I cannot provide sources to the actual data. Ultimately, I need to somehow translate a standard WMS request into an appropriate WFS request. For example, if this is requested by the client:

http://localhost/cgi-bin/mapserv.exe?map=C:/ms4w/apps/test/test.map&SERVICE=WMS&VERSION=1.1.1&LAYERS=WFS_Test&HEIGHT=600&WIDTH=600&REQUEST=GetMap&SRS=EPSG:4326&FORMAT=image/png&BBOX=-104.25,39.0,-104.0,39.25

The below Layer definiton, with the BBOX hard coded as part of the wfs_filter metadata the Mapfile, produces the correct result.  But can I pass the WMS BBOX into the wfs_filter BBOX tag with each request?

    LAYER
      DEBUG 3
      NAME   "DigitalGlobe_WFS"
      TYPE    POLYGON
      STATUS  DEFAULT
      EXTENT  -180 -90 180 90
      PROJECTION
        "init=epsg:4326"
      END
      UNITS DD
      MAXSCALE 144000
      CONNECTIONTYPE WFS
      CONNECTION "https://placeholderservice.com/wfsaccess?
      METADATA
        "wfs_typename"          "WFS_Test"
        "wfs_version"           "1.0.0"
        "wfs_request_method"    "POST"
        "wfs_connectiontimeout" "60"
        "wfs_maxfeatures"       "50"
        "wfs_geometryname"      "geometry"
        "wfs_filter"            "<And>gml:Box><gml:coordinates>-104.25,39.00 -104.00,39.25</gml:coordinates></gml:Box></BBOX><PropertyIsEqualTo><PropertyName>source</PropertyName><Literal>NEW</Literal></PropertyIsEqualTo></And></Filter>
  END


texts: How to create a customizabe table in QGIS Print Composer?

texts:
I would like to add a custom table to a map layout created in QGIS. I would like a fair bit of control over the look of the table (borders, background colours of the cells, etc.), and it appears that the "Add Attribute Table" tool doesn't let you control those things. I tried created what I wanted in Excel, and saving as an image, but of course it's pixelated. I would like a vector-based solution if at all possible. Any thoughts?

I considered creating an .html table, and inserting that; is this the best method?


texts: Getting the actual values from image.normalizedDifference in Google Earth Engine

texts:
I'm trying to display an image after normalizedDifference function in Google Earth Engine
but the default is that negative input values are forced to 0 so that the result is confined to the range (-1, 1).
Can I change the default and get the real values from the calculation?
def addLSWI(image):# add NDVI band to collection
    return image.addBands(image.normalizedDifference([ 'B8','B12',]).rename('LSWI'))

imageCollection = ee.ImageCollection("COPERNICUS/S2");
imageBand_LSWIBand = imageCollection.map(addLSWI)
image_XT = ee.Image(imageBand_LSWIBand.filterDate('2015-10-28','2021-03-30').filterMetadata('MGRS_TILE','equals','36RXT').first());

The results from image_XT.select('LSWI').getInfo():
{'type': 'Image',
 'bands': [{'id': 'LSWI',
   'data_type': {'type': 'PixelType',
    'precision': 'float',
    'min': -1,
    'max': 1},


texts: Add layer to GeoServer from one PostGIS table

texts:
I have a Geonode/Geoserver application with a lot of base layers added to it. I have also developed a Django App plugged into Geonode that lets a user upload shapefiles. These shapefiles are converted using shp2pgsql and stored in a uploads PostGIS table, with each upload having a unique id key.
Once a user uploads a shapefile, the app runs SQL queries to compare (check for overlapping areas) the uploaded shapefile with existing base layers and results are stored in a database table. So far, this has been done.
But now, I need to show the upload shapefile on top of existing base layers in a map viewer (Geonode). Wondering what would be the best approach.
Do I have to add the uploaded shapefiles to GeoServer? Or can I set the uploads table as a PostGIS datastore in GeoServer and have it create layers from each row?


texts: What to do about Geoserver WPS warning" "Could not get a ServiceInfo"?

texts:
This is the warning I am getting when I try to use this WPS service from QGIS,
and the processing time is taking too long, but doesn't terminate.

What might be the problem?

17 Oct 12:33:25 INFO [geoserver.wps] -
Request: getCapabilities
    service = WPS
    baseUrl = http:// localhost:8080/geoserver/
17 Oct 12:33:31 WARN [geoserver.ows] - Could not get a ServiceInfo for service wps
                                       thus could not check if the service is enabled
17 Oct 12:33:31 INFO [geoserver.wps] -
Request: describeProcess
    service = WPS
    version = 1.0.0
    baseUrl = http:// localhost:8080/geoserver/
    identifier[0]:
        value = gt:VectorToRaster


texts: How determine data type of points in geotiff DTM?

texts:
I have some DTMs in geotiff files. I need to read them in C# and I've decided to use libtiff for that purpose. After some investigations, I've found that there are tiled tiffs, so I use ReadTile method to get raw points values. The strange thing for me, newbie to GIS technologies, is that one file seems to storage point elevation values as integers, but another one - as floats (4-byte, single precision). That's the result of comparison of values got by my code and  or  values from *.tif.aux.xml file.
So, my question is - is there any method to determine whether geotiff file stores values as intergers either floats? There is no answer for that in geotiff specs.


texts: ArcGIS Server edit role fails with "Another administrative operation is currently accessing the store"

texts:
I just installed a vanilla instance of ArcGIS Server 10.3.1 on Windows 2008.  From the manager webapp I opened Security -> Settings, then configured the server to use:

User Store: LDAP
Role Store: LDAP
Authentication Tier:    Web
Authentication Mode:    Web

Then I opened Security -> Roles and clicked the edit icon next to the role I want to use as administrator.  I get the expected message:

  Role information is being managed externally in your enterprise's LDAP server. As a result, this information cannot be modified in Manager. Only the 'Role Type' attribute, which is managed by ArcGIS Server, can be modified.

I click the radio button to set Role Type to Administrator and click Save.  After working on it for 30 seconds or so, this message is displayed:

  Another administrative operation is currently accessing the store. Please try again later.

I have tried everything I can think of from trying again later, to a full uninstall/reinstall of ArcGIS server...  Any suggestions as to what could be causing this?

--- UPDATE ---

I tried following this process on a brand new Windows 2012 VM.  Same problem...  The headache continues...

--- UPDATE 2 ---

I uninstalled everything and this time installed 10.4 following the same instructions.  Now when I try to set the Role Type, I get a much more informative error:

    Illegal char <|> at index 52: C:\arcgisserver\config-store\/security/roles\global |-| global arcgis administrators.json\

It looks like ArcGIS LDAP integration does not support all valid names.  I can't really see a way around this other than:

File a report with ESRI
Replicate group to another one with a name not including an | character.

ARGH!

--- Update 3 ---

For ESRI, the role name is Global |-| Global ArcGIS Administrators


texts: How to write a script that re-projects vector datasets in a folder?

texts:
This is probably a novice question, but I am a novice with python.  I am trying to write a script that re-projects vector datasets in a folder.  However there are rules in which to do this.  They are:

Must re-project shapefile vector datasets in the folder to match the
target dataset's projection.
Must append "_projected" to the end of each projected dataset name. For example:
CityBoundaries_projected.shp.
Must skip projecting any datasets that are already in the target projection.
Must report a geoprocessing message of which datasets were projected. In this message, the dataset names can be separated by spaces. In the message, do not include datasets that were skipped because they were already in the target projection.
Must not contain any hard-coded values such as dataset names, path names, or projection names.
Must be made available as a script tool that can be easily run from ArcToolbox by
someone with no knowledge of scripting.

The part I am struggling with is the bolded requirement b/c I cannot hard code the projection I want all of the files to have.  I need it to assign the projection of the remaining shapefiles from that of the targetFC.  So I don't know what to assign to my "outCS" variable. Oh...and I cannot use the batch projection tool.  I have to make a variation of a batch project tool by running the Project tool inside a loop.  Here is a sample of my script. Any help would be greatly appreciated.

import arcpy
import os

#set up the paths
targetFC = r"C:\MGIS\geog485\Lesson2\Lesson2\CountyLines.shp"
folderToExamine = r"C:\MGIS\geog485\Lesson2\Lesson2"

#get spatial reference for the target feature class
targetDescribe = arcpy.Describe(targetFC)
targetSR = targetDescribe.SpatialReference
targetSRName = targetSR.Name

# Get a list of my feature classes
arcpy.env.workspace = folderToExamine
listOfFCs = arcpy.ListFeatureClasses()

#Loop through the list of FCs
for currentFC in listOfFCs:
#Read the spatial reference of the current one
currentFCDescribe = arcpy.Describe(currentFC)
currentFCSR = currentFCDescribe.SpatialReference
currentFCSRName = currentFCSR.Name

if currentFCSRName != targetSRName:
    print "Spatial references don't match"
else:
    print "Spatial references do match"

if currentFCSRName == targetSRName:
    continue
else:
    # Determine the new output feature class path and name
    **outCS =**
    arcpy.Project_management(currentFC, outCS)
    rootName = currentFC[:-4]
    print rootName
    print rootName +"_projected.shp"


texts: How to run Intersect in QGIS when invalid geometry?

texts: I have two layers to make an intersection. I use QGIS 2.18.22 and Vector - Geoprocessing Tools - Intersect - run is not working, because "Input layer A contains invalid geometries. Unable to complete intersection algorithm." Anyone having any idea how to fix it? Layers are here: https://uloz.to/tam/_BGkW3f7iUrDm


texts: Using spatial filters in cross-layer filtering by GeoServer?

texts:
I am working with GeoServer, WFS and have querylayer module installed.

I have been trying to extract information around a featured points. For example, the location of all hospitals within 1000 meters of all schools  over a region. I know I can use DWITHIN on the layer with hospital information are in, and use OR to join multiple DWITHIN
&typeNames=buildings&cql_filter=DWITHIN(geometry,x,y,1000,meters) AND type=hospital

but this requires me to find out the location of the school manually and construct a DWITHIN for every school over the region.

So then I read about cross layer filter here and tried this
&typeNames=buildings&cql_filter=DWITHIN(geometry,collectGeometries(querycollection('buildings','geometry','type=''school'''))) AND type=hospital

This works perfectly if my layer is small, but my layer consist of data for the whole global, so eventually break on querylayer's limits. So I wanted to limited the region for querylayer and did this
&typeNames=buildings&cql_filter=DWITHIN(geometry,collectGeometries(querycollection('buildings','geometry','type=''school'' AND BBOX(geometry,A,B,C,D)'))) AND type=hospital
Then it stopped working. This does not even output errors, just empty dataset.

If spatial filtering can work inside a querylayer, then how would I use it?


texts: RuntimeError: Object: Error in executing tool - ArcPy - ArcGIS Pro - Apply Symbology from Layer

texts:
I have this script to apply the symbology from one layer to many. But when I try and run it, I get the following error: RuntimeError: Object: Error in executing tool

Here is the rest of my code:

import arcpy

# Reference current project and mapframe.
relpath = r'C:\arcGIS_Shared\Python'
p = arcpy.mp.ArcGISProject(relpath + r'\CenterHeatMaps.aprx')
m = p.listMaps('Map')[0]

# Layer file symbology to be imported.
symbology_layer = (r"C:\arcGIS_Shared\Python\Annual2017_CA_ReferenceLyrx.lyrx")

# Layer to receive imported symbology.
layers = m.listLayers("CA201*")
for lyr in layers:
    #Maintain the symbology ranges
#   arcpy.update_symbology('MAINTAIN')
    #Apply the symbology
    arcpy.ApplySymbologyFromLayer_management(layers, symbology_layer)

    #print results
    print(lyr.name)
    count = arcpy.GetMessageCount()
    print (arcpy.GetMessage(count-1))

What am I missing or doing wrong?


texts:
I just updated ArcGIS Desktop (ArcMap) to 10.1 and it is so slow.  When I load my template mxds they can take 30 plus seconds and if I have to change the view it's another 30 seconds.
Does anyone have any ideas or methods to speed up the loading and redrawing purpose?


texts: How to save QGIS Project without loading it with PyQGIS

texts:
I wrote a processing tool and as part of the tool, I am zooming to the selected object(s) and exporting (writing) the current state to a .qgz file. However, it immediately loads this new .qgz  to QGIS and replaces the one which was open when I called the processing tool. I want to save the current state with another name but not load it. Is there any way to do that?
Here is the part of the code where I write .qgz
            # Zoom to Planrahmen
            canvas = iface.mapCanvas()
            initial_canvas = canvas.extent()
            canvas.setExtent(feat.geometry().boundingBox())
            tgtCanvasScale = round(canvas.scale() + 1500, -3)  # rounds up to nearest 1000 and + another 1000
            canvas.zoomScale(tgtCanvasScale)

            # Export to .qgz
            projectName = "GZP_{}_{}_{}_{}_print.qgz".format(gem_name, int(gzpID), labelGZP, karte_bez)
            projectOutTEMP = os.path.join(outDir, projectName)
            project.write(projectOutTEMP)

            feedback.pushInfo("QGZ {} erfolgreich erzeugt.".format(projectName))

            # Zoom back to initial extent
            canvas.setExtent(initial_canvas)
            initCanvasScale = round(canvas.scale())
            canvas.zoomScale(initCanvasScale)


texts: Reordering of layer properties fields in QGIS possible?

texts:
I'm wondering why such simple things are not possible.

I have 5 fields on my layer and the last one shall be on the top now. But there is no way to do this, is this correct?


texts:
I have a great number of area geometries, and line- and point-geometries which have a geometric relation to the area. That means that the distance and the orientation of the lines/points to the areas are important.

Now I got a new version of the areas and found that their position and rotation have changed! The result is, that my lines/points doesn’t fit to the (new) areas. Distances and angles are wrong.

The areas were transformed very individual, so I can’t transform my line/point features with a single transformation. I also can’t use the simple move and rotate functions from QGIS because the line features are often related to more than one area, if I would move it in order to correct the distance to one area I would destroy its geometric relation to other areas.
I can use the really great vertex tool. I select all vertices from the lines which are related to an area …

… and move it from the corner of the old area to the corner of the new area.

In this way I would proceed for each area. But the problem is that the vertex tool only moves my points and vertices. The needed rotation is not possible. You can take the line starting at the corner of the old area as an example. It has a 90 degree angle to the edge of the old area. After transformation it should have a 90 degree angle to the edge of the new area. The selected vertices must be rotated around the corner of the area with an angle (in my case calculated as the difference between the absolute angles of the edges of the new/old areas).

My question is: Anyone have any ideas how I can solve this problem? Is there a plugin for it or do you have some code snippets for me, which brings me on the way to create an own plugin function or can you pimp the vertex tool?


texts: spatial filtering data in R between two vectors

texts:
I have two shapefiles in R. One contains a vector of coastline (called coastline) and the other file of parcels (called parcels). I want to filter the parcel data for only those parcels within 1 mile of the coastline.

How do I find the data only inside the 1 mile coastline radius? Assume it would be something in the rgeos package.


texts: Displaying an icon at each point in list of features using OpenLayers 6

texts:
I am trying to create a layer on a map to indicate a trail in addition to LineString. I want to display an icon each point.
I create my map layer
var  historySource = new ol.source.Vector({
            features: history
        });

        var markersTrailMarkerLayer = new ol.layer.Vector({
            name: 'HistoryTrailIcons' + num, // needed?
            id: 'HistoryTrailIcons' + num,
            source: historySource,
            setZIndex: 500
        });
        map.addLayer(markersTrailMarkerLayer);

I then loop through an array creating a marker for each point:
for (var x = 0 ; x < historyArray.length ; x++)
{
    var lonlat = new ol.proj.transform([historyArray[x].longitude, historyArray[x].latitude], 'EPSG:4326', 'EPSG:27700');

    points.push(lonlat); // for the line string used later

    // Create history marker
    var marker = new ol.Feature({
        geometry: new ol.geom.Point(lonlat),
        id: historyArray[x].num+ "_" + x,
    });
    marker.setId(historyArray[x].num+ "_" + x);

    var iconStyle = new ol.style.Style({
        image: new ol.style.Icon(({
                    anchor: [0.5, 1],
                    src: "/images/historymarker.png",
                    opacity: 1,
                    id: historyArray[x].num+ "_" + x
                }))
            });
            marker.setStyle(iconStyle);

    history.push(marker);
}

It appears to being created but is not visible on top of the line string. Can anyone help?
The LineString I mentioned is a separate layer drawn elsewhere in the code but I have included as requested:
var lineStringFeature = new ol.Feature({
                geometry: new ol.geom.LineString(points),
                name: esn
            });
            lineStringFeature.setId(esn);

            // Add history trail line
            var lineString = new ol.layer.Vector({
                id: "HistoryTrail" + num,
                source: new ol.source.Vector({
                    features: [lineStringFeature]
                }),
                style: new ol.style.Style({
                    stroke : new ol.style.Stroke({
                        color: '#0000ff',
                        width: 5
                    })
                })
            });
            map.addLayer(lineString);


texts: Getting QGIS to use dedicated graphics card (GPU) over integrated?

texts: does anyone know how i'd be able to get QGIS in either windows or Ubuntu to use the dedicated graphics card over the integrated graphics? I am attempting to speed up the processing of a supervised classification algorithm. The algorithm in question  is the semi automatic classification plugin.


texts: missing entry for search statistics geonetwork 3.10.1 / war file installation

texts:
Installed the geonetwork version 3.10.1 from War File (Tomcat 7). Elasticsearch 7.2, Kibana 7.2 running on seperated servers. What we did:
Elasticsearch:
Load the indexes:
curl -O https://raw.githubusercontent.com/geonetwork/core-geonetwork/10.0.1/es/config/features.json

curl -H 'Content-Type: application/json' -X PUT http://host:9200/gn-features -d @features.json

curl -O https://raw.githubusercontent.com/geonetwork/core-geonetwork/10.0.1/es/config/records.json

curl -H 'Content-Type: application/json' -X PUT http://host:9200/gn-records -d @records.json

curl -O https://raw.githubusercontent.com/geonetwork/core-geonetwork/10.0.1/es/config/searchlogs.json

curl -X -H 'Content-Type: application/json' PUT http://host:9200/gn-searchlogs -d @searchlogs.json

Kibana:

Set server.basePath: "/dashboards" #war File is deployed as tomcat ROOT.
Set kibana.index ".dashboards".
Import the export.json from https://github.com/geonetwork/core-geonetwork/blob/master/es/es-dashboards/data/export.json under saved objects in kibana.

GeoNetwork:

Comment <!--<filter-mapping> <filter-name>XFrameOptionsFilter</filter-name> <url-pattern>/*</url-pattern> </filter-mapping>--> in web.xml.
Set parameters in config.properties file for (searchlog,records,features) index, kb.url and es.url.
Set spring.profiles.active value = "es" for dashboard in web.xml.

Health check:

CswGetCapabilitiesHealthCheck: OK

CswGetRecordsHealthCheck: OK
DashboardAppHealthCheck: OK   Dashboard application is running.
DatabaseHealthCheck: OK
DeadlockedThreadsHealthCheck: OK
FieldCacheSanityHealthCheck: OK
Free file handles: OK
FreeConnectionsHealthCheck: OK
HarvestersHealthCheck: OK
IndexHealthCheck: OK   0 records indexed in remote index currently.
LuceneIndexHealthCheck: OK
NoIndexErrorsHealthCheck: OK

but now only content statistics works; the entry for search statistics is missing.
Hope someone here have an idea where my mistake is.
Here the missing search statistics enty:


texts: Remove points that are close to a border when building a grid

texts:
I want to remove the points that are close to a boarder when building a grid, i am using the sf library from R.

library(sf)
library(raster)
library(ggplot2)

shp <- getData('GADM', country = 'aut', level = 0) %>% st_as_sf()

grid <- st_make_grid(shp, n=c(25, 25),
                     what = "centers",
                     square=TRUE) %>% st_intersection(shp)

ggplot() +
  geom_sf(data = shp) +
  geom_sf(data = grid, col ="red")

Example of the problem:


texts: How to map/analyze artifacts that only correspond to 1x1m units?

texts:
The archaeological site has over 4 million lithic (stone tool artifact) remains from 1 occupation level, so individual proveniencing is impossible. Instead, the lithics from each 1x1m unit, were put into a single bag. I have analyzed 100,000 of these lithics from a section of the site, and each one of those lithics have 7 corresponding features. So far everything has been recorded only on paper since I am unsure how to continue.

1) I want to incorporate this lithic analysis data onto an already existing vector based recreation of the site. I will use this to create density maps of different lithic characteristics (by size, type, retouch, completeness, etc.) and if possible do some Exploratory Data Analysis. I am a relatively new student of GIS, and honestly am not sure how to go from my hand written recordings of the lithic data to a spatially referenced and manipulable layer(?).

2) A major part of this project will be to create a workable database of the information transcribed from paper notes. I have 100,000 lithics (with 7 attribute data points each) that are spatially grouped into 1x1m units (~5,000 lithics/unit average). For the actual mapping of this information I understand that I will most likely have to group all related lithics within the same unit together (through pivot tables I think) but is there a way/is there any benefit to individually record all 100,000 lithics as their own feature with its own attribute data, while allowing me still use this information like explained in #1 above?

I know I am asking a lot and some of it is probably very basic, but my university is lacking in GIS.


texts: Determining elevation of polylines using ArcMap

texts:
I am attempting to determine the elevation of about 100k different polylines. I have an associated DEM for my state that is in the same coordinate system.
What tool(s) could I use to determine the mean elevation of each polyline and add this information to the attribute table?


texts: Heads up digitizing of rooms and doors using ArcMap

texts:
I am digitizing the floorplans of various schools in our county. This is for the Fire Department so they want to be able to see where doors are within the building. The doors do not necessarily have to be indicated; a gap in the wall is sufficient.

It would be much more efficient for me to use the 90 degree tool and knock out these classrooms in full.

Is there some sort of eraser feature that might split a line segment in two wherever it intersects?

That way I could basically create a rectangle for the room and erase where the door might be to leave a gap.

I've got a few dozen schools to digitize, and by far the most time consuming process is creating small segments around the doors just to leave that gap.


texts: Make point layer from attributes of polygon layer

texts:
How to make point layer from attributes of polygon layer?

I made two polygon vector layer and have given some attributes to them, after that I labelled them but some of labels from two layers are overlapping.  I tried different ways to label but it is not presentable if I place my labels away from center point. Label is shifting away from center point.

I want to label something like AutoCAD in which if some labels are overlapping I can move that specific label and other labels remains at same location.

So, I thought if I make point layer of labels then I can move them as I wish.


texts: Migrating a Data Directory from 2.1.4 to 2.4.3

texts:
Recently, I got a new server and I installed the latest Geoserver version (2.4.3.).
Now I need to migrate all data (layers, stores, styles, users, workspaces, postgis stores...) from older geoserver 2.1.4 to the new one.
Anyone know how can I do it?
Just move the data directory ?

Best regards


texts: Reclassifying raster analysis (like the one in ArcGIS Spatial Analyst) in QGIS

texts:
I am trying to reclassify raster (DEM) in QGIS. I know how it is done in ArcGIS (Spatial Analyst, reclassify option), but I would like to know how it is done in QGIS?.

The point of my exercise is to set all the values of elevation in DEM in the range from 1200-1600 meters as 1.


texts: QGIS failing to download WMTS capabilities

texts:
Is anyone aware of an issue with WMTS layers in QGIS 3.1.8?

When I create the connection I get the "failed to download capabilities" error. I can add this WMTS into another GIS software just fine and also browsers return the capabilities as expected when requesting
...?request=getcapabilities&service=wmts
I can add WMS from the same server to QGIS with no issues.
It's a local WMS/WMTS so can't share the publicly accessible URL.


texts: Running gdal ogr2ogr command from python

texts:
I have an ogr2ogr command which uploads an ESRI File Geodatabase to a PostgreSQL database, updating/overwriting the existing data:

ogr2ogr -f "PostgreSQL" PG:"host=hosting_url port=5432 dbname=database_name user=username password=password" -lco SCHEMA=schema_name "D:\data\Data.gdb" -overwrite -progress -lco OVERWRITE=YES --config PG_USE_COPY YES

This works as expected when run from OSGeo4W Shell (installed with my QGIS installation). However, I'd like to automate this command to run at the same time as an existing python script I have. So my challenge is to run the ogr2ogr command from python. How can I make it compatible?

I've tried:

import os, subprocess
command = ["C:\\Program Files\\QGIS 2.18\\bin\\ogr2ogr.exe",
        "-f", "PostgreSQL", "PG:\"host=hosting_url port=5432 dbname=database_name user=username password=password\"", "-lco", "SCHEMA=schema_name", "D:\data\Data.gdb", "-overwrite", "-progress", "-lco", "OVERWRITE=YES", "--config", "PG_USE_COPY", "YES"]

I adapted this code from the top answer here: Execute ogr2ogr from python

When I run it I get the error 'returned non-zero exit status 1.'


texts: ArcGIS 10.2.2: How to display attribute values using a representative number of symbols?

texts: I am using ArcGIS 10.2.2 Advanced. I want to produce a map where the value of attributes in each polygon is shown by a representative number of symbols in each polygon. I have used Quantities > Dot Density,  but this assigns the dots to randomised locations within the polygon and I want the symbols to be adjacent to each other as in the example picture attached. Preferably, I would like a solution which allows me to change the dot to a vector symbol of choice as well. I have searched through Stack Exchange and online, but haven't found a way to do this. Has anyone else figured out a way of doing this in ArcGIS?


texts:
When I call the following (using QGIS 3.0.2):

self.layerCatalog =  QgsCore.QgsVectorLayer( "MultiPolygon?%s" % s_fields, nameLayer, "memory")

if self.iface.mapCanvas().isCachingEnabled():
        self.layerCatalog.setCacheImage(None)

I get the error "AttributeError: 'QgsVectorLayer' object has no attribute 'setCacheImage'"

setCacheImage is a function in the parent class QgsMapLayer in QGIS 2. It is not there in QGIS 3. Looking at the Backwards Incompatible Changes site here, I do not see where it was moved or renamed. Anyone know?


texts:
I added a WMTS Layer to a Leaflet map. The tiles are loading but aren't in the correct order.
Has anyone ever encountered this problem? If so, how did you overcome it?

    var map = new L.Map('viewerDiv', {
         crs: L.CRS.EPSG3857
    });
    var url =
        'http://services.arcgisonline.com/arcgis/rest/services/World_Topo_Map/MapServer/WMTS/tile/1.0.0/World_Topo_Map/default/default028mm/{z}/{x}/{y}.png';
    var tilelayer = new L.tileLayer(url, {
        tileSize: 256,
        tms: true,
        attribution: 'World Map'
    });
    map.addLayer(tilelayer);
    map.setView(L.latLng(50, 8), 4);

http://services.arcgisonline.com/arcgis/rest/services/World_Topo_Map/MapServer/WMTS/1.0.0/WMTSCapabilities.xml


texts: How to get layers from a specific store of geoserver using python gsconfig?

texts:
I'm using gsconfig python library to acess geoserver rest, my question is how can i get layers from a specific store ?

I'm tried as follow:

from geoserver.catalog import Catalog

cat=Catalog("http://localhost:8080/geoserver/rest")

layers=cat.get_layers('storename')

but it's giving all published layers from geoserver not from that particular store.

Please help me?


texts: Model builder: select feature based on attribute, iterator?

texts:
This is the data I have:
shapefile with polygons
shapefile with line features

What I eventually want:
line features attributes joined to the polygons that they fall in.
I need to clip each polygon with the lines separately.

These are the steps I want to do:

select polygon from a shapefile based on a name
use this selected polygon to clip several line features
then use spatial join to join these clipped lines to the shapefile with all the polygons based on the location
apply summary statistics to get lines length summed per polygon

I can do this for one polygon..the question I have is how can I tell the model builder that has to do this whole process for each polygon. So I guess I have to insert an iterator in front of the select by attribute tool..but I don't really understand how.
Also, it would be great to get the output in one summary statistic and not in hundred different ones...
is that possible?

thanks so much
sandra


texts: Geoserver SQL View does not pass multiple viewParams value from WMS to PostgreSQL function

texts:
In my application, I am rendering Geoserver layer on map. Base database for Geoserver is PostGresql.

Below give is the request that call WMS layer. You can see, I am passing multiple viewParams as key value(layerid,centerlat and centerlong)

http://localhost:8080/geoserver/wms?viewparams=layerid:9;centerlat:39.44043541908485;centerlong:-83.66638183593751&service=WMS&request=GetMap&layers=myWorkSpace%3Avw_GetLayerData&styles=&format=image%2Fpng&transparent=true&version=1.1.1&tiled=true&height=256&width=256&srs=EPSG%3A3857&bbox=-9549125.069610499,4618019.500877211,-9470853.552646479,4696291.017841231

In Geoserver, I have created a view that run the below given script

select geom from fn_getlayerdata('%layerid%','%centerlat%','%centerlong%')

fn_getlayerdata (integer,float,float): This is a function created in PostgreSQL.

Problem is that fn_getlayerdata function is not called. But If change the Geoserver view as given below it works.

select geom from fn_getlayerdata(%layerid%,88.78,-90.89)

function fn_getlayerdata returns table. Strange thing is that in fn_getlayerdata function value of centerlat and centerlong is not stored in tempcenterlocation table. Also there is no exception.

   CREATE OR REPLACE FUNCTION fn_getlayerdata(
                  layerid integer,centerlat float,centerlong float)
              RETURNS TABLE(geom geometry)
              LANGUAGE 'plpgsql'
           VOLATILE
        AS $BODY$
       DECLARE exec_error VARCHAR;
      BEGIN

INSERT INTO tempcenterlocation(centerlat,centerlong)
values(centerlat,centerlong);
RETURN QUERY
select geom from myTable where dataid>=layerid;

 EXCEPTION WHEN Others THEN exec_error:=SQLERRM;
 Insert Into public."logtable" ("id","errormessage","idtype") VALUES
      (10,exec_error,'fn_getlayerdata');
 END;

    $BODY$;

       ALTER FUNCTION fn_getlayerdata(integer,float,float)
           OWNER TO postgres;

Any Idea What wrong I am doing?


texts:
I have a line shapefile of roads and a point shapefile of car collisions, and am trying to count the number of collisions per road. The points do not exactly intersect with the lines. I'm hoping to end up with a new column in the "roads" attribute table with the number of collisions. I've found a post on how to do that in ArcGIS here: How to produce a 'count' of points intersecting lines?, but haven't been able to figure it out in QGIS.

I've found suggestions to create a buffer of the roads in order to create a polygon and use the Points in Polygon analysis, but the roads data is quite large and creating a buffer was taking ages for me, so I'm hoping there's another option. I've tried using the "Join attributes by location" function, with "roads" as the target layer and "collisions" as the join layer, but I don't know how to set the Precision input so that all of the collision points are counted.  In Arc "Join based on spatial location" has "closest" option when joining points to lines, and I'm hoping that QGIS offers something comparable.

I'm using QGIS 2.18.2