texts: Heatmap of a large polyline dataset in QGIS

texts:
To summarize: I'm having difficulties with creating a heatmap from a large (90021-item) polyline dataset in QGIS 3.10.1.

First, it is worth noting that this question is somewhat similar to Doing heatmap from line data in QGIS? and to Create line heatmap with QGIS. In particular, I have a single polyline (vector) layer containing all 90021 features, each with a few thousand vertices. These lines represent Google Maps-style travel routes, and each one also has an attribute that I'd like to use as a weighting factor in my heatmap.

I initially attempted to follow roughly the procedure described in the third answer to the first question above, which consists of:

Densifying each line (using the native QGIS algorithm for this), so that vertices occur at regular intervals;
Extracting these vertices (again, with the native QGIS algorithm), isolating them from their lines while preserving the weighting attribute; and
Editing the symbology for this vertex layer so that it renders as a heatmap, using the desired attribute as a weight factor.

Unfortunately, the size of my dataset makes this approach infeasible - an earlier dataset, at half the size of my current one, quickly filled my computer's RAM and froze it during this process.

I also tried the approach outlined in the accepted answer to the second question above, using feature blending modes to attain the desired effect. Again, the size of my dataset was problematic. The Addition blending mode, seems to be the best option for heatmap-like output, as it simply adds the RGB components of the pixels to be overlaid. (More layered lines = a lighter color.) However, so many lines overlap that even with the feature color set to 1, 1, 1 (in RGB form), a large portion of the layer becomes white - in effect, there are so many lines overlaid on each other that the RGB values hit their maximum too soon.

Does anyone know how I might be able to proceed? (If it helps, here's the list of QGIS blending modes.)


texts: Finding nearest line/point from polygon using PostGIS?

texts:
I have got the point and the polygon (or line). I have to find the closest coordinate (it can be the "virtual" point on edge) from polygon due to the point.

Do you know how to find out this coordinate?


texts: gdal_merge.py not properly assigning NoData values

texts:
I am trying to create a 4-band raster image using gdal_merge.py but only the first band is getting the correct no data values assigned, according to ArcGIS and PCI. Here is my command:
gdal_merge.py -co COMPRESS=LZW -co BIGTIFF=YES -ul_lr %s -ot Int16 -separate -n -999 -a_nodata -999 -o %s %s' % (extent, outputFile, inputFilesString)

where:
inputFilesString = "input1.tif input2.tif input3.tif input4.tif"

Where there is NoData, ArcGIS is showing NoData for Band1, but 0's for Bands2-4. However, this problem doesn't occur in QGIS (where all bands show NoData where they should). Also, according to gdalinfo (and ArcGIS/QGIS properties), each band's NoData value is correctly -999. All these things were leading me to think it could just be a bug in ArcGIS/PCI.
However, when I read these bands (2-4) into individual arrays and write each of them each back to their own single band .tif (not setting the ND val) using GDAL API, it is writing these NoData areas as 0 while for band1 it correctly writes them as -999. So it is likely an issue in GDAL, not Arc/PCI.
I double checked (in ArcGIS and QGIS), and each of the individual tiffs that are the inputs for gdal_merge have NoData values of -999 that are correctly in areas where there should be NoData.
This is just not translating to the tiff output from gdal_merge.
Is this a bug with GDAL? Or am I not correctly assigning the NoData values in the gdal_merge.py command? From the documentation, it appears the NoData values should be assigned once for all bands as opposed to explicitly defining the values for each band, but I could be wrong.
UPDATE: I reran the exact same command as above but without the "-n -999" string, as suggested in an answer below, but nothing changed.


texts: Leaflet : Change marker icon from onEachFeature instead of pointToLayer

texts:
I'm referencing external data, matching it up with particular fields in feature.properties and using that to decide which layer to put the feature in, all of this from inside the onEachFeature part.

onEachFeature: function(feature, layer) {
...
if (condition1)
  layer.addTo(layer1);
else if (condition2)
  layer.addTo(layer2);
else ...
}

I want to set different icons / markers also in the same process. But I see that custom markers are usually set in the pointToLayer part, like:

pointToLayer: function (feature, latlng) {
  ..
  if(condition1) return L.marker(latlng, {icon: icon1});
  else if(condition2) return L.marker(latlng, {icon: icon2});
  else ...
}

(assume that icon1 and icon2 have already been defined as new L.icon({ .... }) )

It's going to get really repetitive and processor/network-consuming to reference the external data all over again in pointToLayer. Is there a way I can define the feature's icon from the onEachFeature part itself and not have to go through pointToLayer ?

Edit: And even the other way around : Is it possible to do the things you'd normally do in onEachFeature, in pointToLayer instead? The layer argument isn't there in most uses.. can it be passed in? Does leaflet loop through each feature for pointToLayer the same way it does for onEachFeature ?


texts:
Trying to add some Transparency / Opacity to my GeoTiff raster files for a better rendering option when overlapping.

Without transparency

With transparency 20%

Is there any gdal command to achieve this?

Driver: GTiff/GeoTIFF
Files: Gb_127_d_D_0006_georefYCBCR_CJPEG.TIF
       Gb_127_d_D_0006_georefYCBCR_CJPEG.tfw
Size is 14013, 13566
Coordinate System is:
PROJCS["CH1903 / LV03",
    GEOGCS["CH1903",
        DATUM["CH1903",
            SPHEROID["Bessel 1841",6377397.155,299.1528128000008,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[674.4,15.1,405.3,0,0,0,0],
            AUTHORITY["EPSG","6149"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4149"]],
    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
    PARAMETER["latitude_of_center",46.95240555555556],
    PARAMETER["longitude_of_center",7.439583333333333],
    PARAMETER["azimuth",90],
    PARAMETER["rectified_grid_angle",90],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",200000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","21781"]]
Origin = (541546.758578460663557,152626.253885675920174)
Pixel Size = (0.031810463058218,-0.031810463058218)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  541546.759,  152626.254) (  6d40'39.68"E, 46d31'25.15"N)
Lower Left  (  541546.759,  152194.713) (  6d40'39.88"E, 46d31'11.17"N)
Upper Right (  541992.519,  152626.254) (  6d41' 0.60"E, 46d31'25.29"N)
Lower Right (  541992.519,  152194.713) (  6d41' 0.79"E, 46d31'11.31"N)
Center      (  541769.639,  152410.484) (  6d40'50.24"E, 46d31'18.23"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 7007x6783, 3504x3392, 1752x1696, 876x848
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 7007x6783, 3504x3392, 1752x1696, 876x848
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 7007x6783, 3504x3392, 1752x1696, 876x848
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 7007x6783, 3504x3392, 1752x1696, 876x848
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 7007x6783, 3504x3392, 1752x1696, 876x848
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 7007x6783, 3504x3392, 1752x1696, 876x848
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 7007x6783, 3504x3392, 1752x1696, 876x848

EDIT : Proposed workflow to enable transparency on raster data containing no_data values

In my case I had RGBA rasters with no-data values:

Explode bands into one rgb file and one transparency file:

os.system("gdal_translate -b 1 -b 2 -b 3 " + infile + " " +  rgbfile)
os.system("gdal_translate -b 4 " + infile + " " +  transfile)

Compute desired transparency

os.system("gdal_calc -A " + transfile + " --outfile="+ transfilemoded + " --calc=(A/5*4)")
# A was 255 in my case for all data values, 0 for no_data values
# Hint: since we're in Byte, 255*5/4 = 50 and 255/5*4 = 204

Merge the produced rgb and new transparency file

os.system("gdal_merge " + rgbfile + " " + transfilemoded + " -o " + rgbafile + " -separate")


texts: Hydrology Analysis via SAGA - what is the unit?

texts:
I did a couple of analysis via SAGA about an area in Iceland...I've never worked with SAGA before, could anyone tell me  what Units the numbers in the legends for "Catchment Area Recursive" and "topographic Wetness index" are?

Plus: why is the Strahler legend with decimals? shouldn't that be natural numbers?


texts: Scripting (PyQGIS) the match to symbols from file (XML) in QGIS

texts: In QGIS (2.1x) > layer properties > having symbology by categories and symbols classes configured in a XML file (external source) - being "ID_UCF" the key field for the classification - is it possible to automate the operation/option of match to symbols from file? We are investigating around PyQt to try a solution... Any help?


texts: Memory requirements to use TileMill efficently

texts:
When trying to run TileMill with OSM Bright all I get is a "code 3" error, which apparently stands for "out of memory" error.

I'm trying to work on France, but I could restrict the work on one region at a time if it helped with memory.

The PostGIS db runs in a VirtualBox machine running Debian.

My PC is an Intel i5 with 4GB of RAM.

Should I:

buy a new machine (I'd rather not)?
rent something on AWS?
run in a non-interactive way?

Or is there a commercial alternative to TileMill that's more efficient? Any other idea?


texts: Calculate distance between points and nearest polygon in R

texts:
I have a data frame with among others longitude and latitude as variables for about 17,000 points. I also have a shapefile consisting of multiple polygons. The points fall inside or outside the polygons. I would like to add to the data frame the distance to nearest polygon for each point but i have been struggling quite a lot.

This is what i have so far after reading multiple posts on the topic but this only works for points inside the polygons apparently.

library(spatstat)
library(maptools)
shp_Poly <- readShapeSpatial("POLYGON.shp")
W <- as.owin(shp_Poly_BAT)
W <- as.psp(W)
p <- ppp(lon, lat, W)
AVG <- nncross(p, W)

How can i improve this?


texts: Can't interact with map inside a control

texts:
My application has a main window and other controls which I open inside the window, one of these controls is one which contains a map, it appears but I can't move it by holding a dragging also the click event on the map is not provoked.
What is causing this?

XAML:

<UserControl x:Class="ModernUINavigationApp1.Pages.KasrOkeelPages.Inputs"
         xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
         xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
         xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006"
         xmlns:d="http://schemas.microsoft.com/expression/blend/2008"
         xmlns:mui="http://firstfloorsoftware.com/ModernUI"
         xmlns:esri="http://schemas.esri.com/arcgis/client/2009"

         mc:Ignorable="d"
         d:DesignHeight="300" d:DesignWidth="300" MouseDoubleClick="UserControl_MouseDoubleClick" MouseClick="MyMap_MouseClick" Focusable="True">
<Grid Style="{StaticResource ContentRoot}" x:Name="LayoutRoot">
    <esri:Map x:Name="MyMap" UseAcceleratedDisplay="True" Cursor="Cross" Extent="330671.7091,2854689.5845,344178.0091,2862125.3845"  IsEnabled="True" DragOver="MyMap_DragOver" IsLogoVisible="False" WrapAround="True">
        <!--Local Dynamic Layer-->
        <esri:ArcGISLocalDynamicMapServiceLayer ID="MYLAYERS" Path="C:\Users\Ahmad Elsayed\Desktop\test app\GIS APP2\packTest\Untitled2.mpk"/>

    </esri:Map>
</Grid>

C# :

    using FirstFloor.ModernUI.Windows.Controls;
    using System;
    using System.Collections.Generic;
    using System.ComponentModel;
    using System.Linq;
    using System.Text;
    using System.Threading.Tasks;
    using System.Windows;
    using System.Windows.Controls;
    using System.Windows.Data;
    using System.Windows.Documents;
    using System.Windows.Input;
    using System.Windows.Media;
    using System.Windows.Media.Imaging;
    using System.Windows.Navigation;
    using System.Windows.Shapes;

    namespace ModernUINavigationApp1.Pages.KasrOkeelPages
    {
        /// <summary>
        /// Interaction logic for Inputs.xaml
        /// </summary>
        public partial class Inputs : UserControl
        {
            public Inputs()
            {
                InitializeComponent();
            }

            private void MyMap_DragOver(object sender, DragEventArgs e)
            {

            }

            private void UserControl_MouseDoubleClick(object sender, MouseButtonEventArgs e)
            {

            }

private void MyMap_MouseClick(object sender, ESRI.ArcGIS.Client.Map.MouseEventArgs e)
        {

        }
        }
    }


texts: 3D OSM buildings seems to be floating and shifted after adding terrain effect in Cesium

texts:
I had made a 3D view of the osm buildings using cesium. In which i had added a terrain effect using the standard cesium terrain provider as below

var viewer = new Cesium.Viewer('cesiumContainer');
var terrainProvider = new Cesium.CesiumTerrainProvider({
  url : '//assets.agi.com/stk-terrain/world'
});
viewer.terrainProvider = terrainProvider;
var dataSource2 = new Cesium.GeoJsonDataSource();
var promise = dataSource2.load('../../SampleData/Hubli_Buildings.geojson');
promise.then(function(dataSource2) {
    viewer.dataSources.add(dataSource2);
    viewer.zoomTo(dataSource2);

    //Get the array of entities
    var entities2 = dataSource2.entities.values;

    var colorHash = {};
    for (var i = 0; i < entities2.length; i++) {
        //For each entity, create a random color based on the state name.
        //Some states have multiple entities, so we store the color in a
        //hash so that we use the same color for the entire state.
        var entity = entities2[i];
        var name = entity.Elev_in_m;
        var color = colorHash[name];
        if (!color) {
            color = Cesium.Color.BROWN;
            colorHash[name] = color;
        }

        //Set the polygon material to our random color.
        entity.polygon.material = color;
        //Remove the outlines.
        entity.polygon.outline = false;

        //Extrude the polygon based on the state's population.  Each entity
        //stores the properties for the GeoJSON feature it was created from
        //Since the population is a huge number, we divide by 50.
        entity.polygon.extrudedHeight = entity.properties.Elev_in_m;
        entity.polygon.height = entity.properties.Elev_in_m;
    }
}).otherwise(function(error){
    //Display any errrors encountered while loading.
    window.alert(error);
});

After adding the terrain, my geojson 3D seems to be floated and shifted as below

I would like to know how to avoid this shifting effect after adding the terrain.


texts: Creating a Grid in QGIS of with 1meter intervals

texts:
I need to create a grid over my QGIS shapefile, with 1m intervals.
I know the procedure of creating the grid; view -> decorations -> grid

My question is, what is the x and y intervals that I need to provide in order to get a grid with 1 meter intervals, when I am creating a line grid?


texts: Converting files in a specific folder specified by a relative path

texts:
I am writing a script to automate the process of converting a group of rasters, the fact is that I need to set a relative path, it will be used on different pc. But there is always an error [error-000732], how to solve this problem?
The code that I have now
error stack:
Traceback (most recent call last):
  File "C:\projects\prototype\Конвертация.py", line 15, in <module>
    arcpy.RasterToOtherFormat_conversion("input_files", "output", "TIFF")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\conversion.py", line 2658, in RasterToOtherFormat
    raise e
ExecuteError: РќРµ СѓРґР°РµС‚СЃСЏ РІС‹РїРѕР»РЅРёС‚СЊ. РџР°СЂР°РјРµС‚СЂС‹ РЅРµ РєРѕСЂСЂРµРєС‚РЅС‹.
ERROR 000732: Input Rasters: РќР°Р±РѕСЂ РґР°РЅРЅС‹С… input_files РЅРµ СЃСѓС‰РµСЃС‚РІСѓРµС‚ РёР»Рё РЅРµ РїРѕРґРґРµСЂР¶РёРІР°РµС‚СЃСЏ
РќРµ СѓРґР°Р»РѕСЃСЊ РІС‹РїРѕР»РЅРёС‚СЊ (RasterToOtherFormat).

import os
import glob
import arcpy
from arcpy import env
inpath = r".\\data\\ice\\input"
env.workspace = os.path.realpath(inpath)
data_rasters = arcpy.ListRasters("", "BIL")
for raster in data_rasters:
    desc = arcpy.Describe(raster)
    path = desc.path
    input_files = str(path) + "/" + raster
    arcpy.RasterToOtherFormat_conversion("input_files", "output", "TIFF")


texts: Problem with autofilling output layer names while batch processing

texts:
I've been using QGIS "Raster boolean OR" to batch process a set of rasters (using LTS version 3.16.6). This has been working as hoped/expected, with the exception that I've not been able to get autofill to use the input layer names when generating the output layer names. When I select the raster layer inputs, they populate the "Input layers" column as a list of "1 inputs selected"-elements, which doesn't seem to provide a readable parameter when using the autofill settings:

When trying to use the autofill function to generate the output layer names using the input layer names it throws me a python error (which I think is saying it's unable to read the list of inputs):

TypeError: unhashable type: 'list'
Traceback (most recent call last):
File
"C:/PROGRA~1/QGIS3~1.16/apps/qgis-ltr/./python/plugins\processing\gui\BatchOutputSelectionPanel.py",
line 112, in showSelectionDialog
if v in QgsProject.instance().mapLayers(): TypeError: unhashable type: 'list'

I've been unable to find any clues as to how to work around this problem. I haven't even seen any examples of people experiencing similar issues.
Anyone have any ideas?


texts: Which fields are standardized to describe a ShapeFiles feature?

texts:
I would like to store data into a database which is available here either as JSON or as ShapeFiles. The JSON looks something like this:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "id": "HUNDEZONEOGD.12884",
      "geometry": {
        "type": "Point",
        "coordinates": [
          16.338031572401313,
          48.20739110196646
        ]
      },
      "geometry_name": "SHAPE",
      "properties": {
        "PARK": "Lerchenfelder Gürtel 30-32",
        "FLAECHE": "463 m²",
        "TYP": "Hundezone",
        "TELEFON": "01 4000 8042",
        "EINFRIEDUNG": "Stabilgitterzaun 1m",
        "WEBLINK1": "http://www.wien.gv.at/umwelt/parks/hundezonen.html"
      }
    },
    ...
  }
],
"crs": {
  "type": "EPSG",
  "properties": {
    "code": "4326"
  }
},
"bbox": [
  16.26757004155759,
  48.13103267739301,
  16.480528967294525,
  48.28714492813137
  ]
}

I also downloaded the ShapeFiles alternative and tested different Ruby libraries to extract the information.

zhm/gdal-ruby
dazuma/rgeo with dazuma/rgeo-shapefile
nofxx/georuby

All of them are capable of reading the properties of the actual feature. I noticed however that they do not support reading the following fields:

feature.id
feature.type
feature.geometry_name
feature.geometry.type // supported by RGeo

Again I looked into the data using Quantum GIS and realized that the ShapeFiles do not contain the mentioned fields (feature.id, feature.type, feature.geometry_name). Does this mean these fields are not part of the ShapeFiles standard? I am confused why the publisher includes the fields in the JSON but not in the ShapeFiles.

I found out how each of the libraries accesses the attribute feature.geometry.type:

zhm/gdal-ruby: feature.get_geometry_ref.get_geometry_name
dazuma/rgeo: feature.geometry.geometry_type.to_s
nofxx/georuby feature.geometry.text_geometry_type

It was a bit hard to find with gdal-ruby and rgeo.


texts: How to apply svm to a dataset with a numerical (not categorical) dependent variable?

texts:
Context: I have some band values from a sentinel - 2 derived .tiff I now want to make a prediction regarding areas where i have no actual field data i.e. make a carbon map of sorts

Libraries are: library("e1071") library("caret")

IndependentVar<-AGCplot[,3:19] DependentVar<-AGCplot[,2] mod_svm <- svm(AGC~B2+B3+B4+B5+B6+B7+B8+B8A+B11+B12+NDVI+IRECI+VV+HV+NDWI+MMRI+MASI, data=AGCplot) summary(mod_svm) pred <- predict(mod_svm,IndependentVar)

I can't run a confusion matrix as ways of evaluating the model because the dependent variable is numerical. What can i run to valuate how my model did?


texts:
I can't find what is wrong.  I try write an style on GeoServer.

This is not ok.:

    <FeatureTypeStyle>
                <FeatureTypeName>Odcinki</FeatureTypeName>
          <Rule>
      <Name>Magistrala A</Name>
        <Title>&lt; 1700</Title>
  <MinScaleDenominator>1700</MinScaleDenominator>
  <MaxScaleDenominator>6771</MaxScaleDenominator>
     <ogc:Filter>**line.24**
        <ogc:PropertyIsEqualTo>
         <ogc:PropertyName>mpec_magis</ogc:PropertyName>
         <ogc:Literal>1</ogc:Literal>
        </ogc:PropertyIsEqualTo>
    </ogc:Filter>
        <LineSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#009933</CssParameter>
          <CssParameter name="stroke-width">4</CssParameter>
        </Stroke>
        </LineSymbolizer>
  </Rule>
  <Rule>  <!-- thick line drawn first-->
    <MaxScaleDenominator>6771</MaxScaleDenominator>
            <ogc:Filter>**line.40**
        <ogc:PropertyIsEqualTo>
         <ogc:PropertyName>mpec_magis</ogc:PropertyName>
         <ogc:Literal>1</ogc:Literal>
        </ogc:PropertyIsEqualTo>
    </ogc:Filter>
    <LineSymbolizer>
      <Stroke>
          <CssParameter name="stroke">#009933</CssParameter>
          <CssParameter name="stroke-width">2</CssParameter>
      </Stroke>
    </LineSymbolizer>
        </Rule>

I have got errors:

line 24: cvc-complex-type.2.4.a: Invalid content was found starting with element 'ogc:Filter'. One of '{"http://www.opengis.net/sld":Symbolizer}' is expected.

line 40: cvc-complex-type.2.4.a: Invalid content was found starting with element 'ogc:Filter'. One of '{"http://www.opengis.net/sld":Symbolizer}' is expected.

I don't understand what is wrong, beyond it is LineSymbolizer.

My previous wersion without ScaleDenominator is ok.:

    <FeatureTypeStyle>
 <Rule>
   <Name>Magistrala A</Name>
   <ogc:Filter>
     <ogc:PropertyIsEqualTo>
       <ogc:PropertyName>mpec_magis</ogc:PropertyName>
       <ogc:Literal>1</ogc:Literal>
     </ogc:PropertyIsEqualTo>
   </ogc:Filter>
   <LineSymbolizer>
     <Stroke>
       <CssParameter name="stroke">#009933</CssParameter>
       <CssParameter name="stroke-width">2</CssParameter>
     </Stroke>
   </LineSymbolizer>
 </Rule>

For any suggestions many thanks.
I greet


texts:
I have 111 Landsat8 scenes which I want to calculate several indices for.
The algebra to calculate them is most of the time Index= (b1-b2)/(b1+b2).
I have the single bands (bands 1-11 respectively) all stored in one folder, named LC08_L1TP_194023_20140225_20170425_01_T1_B1 etc.

I am looking for a good way to automatically calculate the indices.

I have made an approach in R already, which works fine for one file, but I don't know how to make a loop reading the correct files in the folder out of it.

red=raster("Path:/LC08_L1TP_194024_20171015_20171024_01_T1_B4.tif")
nearinfrared=raster("Path:/LC08_L1TP_194024_20171015_20171024_01_T1_B5.tif")

ndvi = (nearinfrared - red) / (nearinfrared + red)
writeRaster(ndvi, filename="Path/NDVI_Test",format="GTiff")


texts: Error in opening rasterfile with ArcView 3.2

texts:
I'm working the first time with ArcView 3.2, but have worked with ArcMap Desktop. I did all the preprocessing of my data in ArcMap, but need a tool (AVErosion) that is just available for ArcView.
My problem is, that all the rasterfiles (.tif) I created just appear as a black blank square in ArcView. In ArcMap and other programs they appear just fine and with information in it.
Does anyone have a clue? I already got the extension 'TIFF 6.0 Image Support' activated.


texts: How to split single lines to represent two-way streets?

texts:
I am working on a project using city road data. The data has every road in the city visualized as lines and with directional data stored in the attribute table (one-way to-from, one-way from-to, and all). I'm wondering if there is a way to split the two-way lines into two separate one-way lines offset from one another but ending in the same point. Any help would be tremendously appreciated!

Thank you!


texts:
Is there any work around to mitigate the amount of time QGIS takes when loading projects that have large tables from postgis layers .
I have tried to turn off this large layers so they are not in use , save , close project and then reopen project hoping that the turned off layers would not load on project start again. But nothing it seems that it will always load the tables in memory even though they are turned off.
Any suggestion ???
using 2.8.2 - Wien on windows


texts: Advice on purchasing new scanner / plotter?

texts:
So I have been tasked with finding options to replace my GIS department's recently deceased HP Designjet 815mfp (all-in-one scanner/plotter/touch screen).  It was used mostly to scan/print plats, and some older maps from archives too.

We would like similar capabilities with a new solution; fast scan/printing, supports up to 42" media, easy to use.  Affordable consumables like ink/printheads would be nice(one can dream). I have been doing research, but there is a bewildering amount of companies and products and minimal reviews.  Budget under $20,000 US.

Are you fellow GIS pros still using HP?  They are so expensive, anyone using other companies, Canon, Epson, Contex?  Would like to hear you recent purchases and experiences Thanks!


texts:
I use this comand

  lpad( aggregate('FIRELAYER','count',"CITYwhithFIRE",(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('0 day')))
+(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('1 day')))
+(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('2 day')))
+(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('3 day')))
+(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('4 day')))
+(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('5 day')))
+(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('6 day')))
+(to_date("DateHourWhenfirehappens")=to_date(to_date(@atlas_pagename)-to_interval('7 day')))),4,0)

for call last seven days of fire points on QIGS Atlas, I would like to know, if there is another simpler way to do this same command.

Table example

The resul is a number of fire, on last seven days
from the day of atlas page name


texts: numpy.unique: ValueError: all the input arrays must have same number of dimensions

texts:
I'm trying to get the unique values of an array but am getting an error that doesn't make any sense.

Here is my pseudo code:

import numpy as np
import gdal

# I have a .tif file, call it inputRaster:

ds = gdal.Open(inputRaster)
dt = ds.GetRasterBand(1).DataType
ncols = ds.RasterXSize
nrows = ds.RasterYSize
outArr = ds.GetRasterBand(1).ReadAsArray()

print dt # prints 2 (uint16)
print ncols # prints 55500
print nrows # prints 55500

print outArr.shape # prints (55500, 55500)
print outArr.dtype # prints uint16

but:

print np.unique(outArr)

gives me the following error:

File "/usr/lib/python2.7/dist-packages/numpy/lib/arraysetops.py", line 197, in unique
 flag = np.concatenate(([True], ar[1:] != ar[:-1]))
ValueError: all the input arrays must have same number of dimensions

My only input is the one array that I want to get unique values for. Why am I getting this error?


texts: Calculating statistics for an Image Service

texts: I haven't been able to find much information at all in regards to ArcGIS 10.1 Function Template Chains... I need to be able to load a .rft in my image service and calculate statistics... so I am using the Arithmetic Function, multiplying the raster by 3.28084 in order to convert from meters to ft. Now this produces the right values, however it does not calculate statistics in order to display the image within the proper range... The pixels calculate, but the image remains symbolized in the original range.  I can manually calculate statistics in ArcCatalog, but if I am supposed to be able to export this function as a .rft for use in my image service, how do I also make it calculate statistics so the user can apply the function and view the image properly??  It seems as if there should be another function called Update Raster.. or Calculate Statistics so when you apply the function it automatically does it, but no such function exists...


texts:
I'm relatively new to using QGIS and GIS platforms in general.

I'm following a guided exercise in QGIS and need to use the 'Spatial Query Tool', however I can't seem to find it anywhere in QGIS 3.4.


texts: Performing statistical analysis of buffers?

texts: I am doing a buffer analysis of incinerators and looking at the output areas within those buffers, is there a way of doing a statistical analysis of several buffers to see if there is a significant pattern between the data within the buffers?


texts: Generating surfaces of different heights from LIDAR points

texts:
I am seeking a way to generate surfaces (polygons in a 3D space) of different heights from LIDAR point clouds, whose spacing is 4 points per square meter. The typical application is like roofs, which may include holes, due to different heights.

However, most tools, which I can find, only generate the boundaries or footprints.

ArcGIS: the examples are like: http://desktop.arcgis.com/en/arcmap/10.3/manage-data/las-dataset/lidar-solutions-data-area-delineation-from-lidar-points.htm and Creating boundary polygon shapefile from set of LAS files using ArcGIS Desktop?
LAStools lasboundary.exe at http://www.cs.unc.edu/~isenburg/lastools/
Python mayavi package: the example is like: https://stackoverflow.com/questions/33376657/from-point-cloud-to-surface-using-python

My initial thought is to do conversion from grid to polygons. Las file can be easy to convert to TIN, and then to grid, but grids are missing height information. Thus, the way is difficult to delineate detailed roof structures.

I also have a 30-centimeter and a 1-meter resolution orthophotos, if these do help, but no stereo photos.

The concrete expectation is like the following picture. The multi-color points are LIDAR point cloud, and the pink lines, which are 'digitized manually', indicate our expectation. However, it is easy to get the outside building boundary (the selected blue-green line), once buildings are classified, but no roofs with different heights (smaller inner polygons) are delineated.


texts: Export collection to individual images in GEE

texts:
I want to clip a MODIS image collection to an area defined by a bounding box in Google Earth Engine and then export the images separately to my google drive. I want to use a map function to avoid for loop, but the export module parameters can not be dynamically set. For example, for the description (that would be the file name), I want it to be equal to the id of each individual image. However, export module expect a string for name and can not accept a computed object. getInfo() also doesn't work in mapped functions.

var imageCol = MODIS.select('LST_Day_1km').filterDate(startDate, endDate);
var clippedCol = imageCol.map(function(image){
    var clipped = image.clip(Bbox);
    Export.image.toAsset({
       image: clipped,
       description: image.metadata('system:index'),
       region: Bbox,
       scale: 1000
    });
    return clipped   });

Any clue?


texts: Loading ADF file from National Elevation Dataset into GRASS with Python?

texts:
I'm trying to load data from the National Elevation Dataset into GRASS/Python, using this Python script (outside of GRASS):

# Imports
import os
import grass_session
from grass.script import setup as gsetup
import grass.script as gscript
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import imagery as i
from grass.pygrass.modules.shortcuts import display as d

# Settings
gisdb = "/home/user/research/grassdata"
location = "USA"
mapset = "PERMANENT"

# Key directories
project = "/home/user/research/radio/analysis/gis"

# Initialization
gisbase = os.environ['GISBASE']
gsetup.init(gisbase, gisdb, location, mapset)

# Sample GIS data
adf_file = os.path.join(project, "sample", "grdn48w095_13", "w001001.adf")
output_file = os.path.join(project, "output", "sample")

print("Importing raster data... ", end = "")
r.in_gdal(input = adf_file, output = output_file, flags = "e", overwrite = False)
print("Done.")

But I get the following error:

Traceback (most recent call last):

  File "/home/user/apps/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2961, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)

  File "<ipython-input-7-77fc8fd256bb>", line 31, in <module>
    r.in_gdal(input = adf_file, output = output_file, flags = "e", overwrite = False)

  File "/usr/lib/grass72/etc/python/grass/pygrass/modules/shortcuts.py", line 47, in __getattr__
    return self.cls('%s.%s' % (self.prefix, name.replace('_', '.')))

  File "/usr/lib/grass72/etc/python/grass/pygrass/modules/interface/module.py", line 475, in __init__
    tree = fromstring(self.xml)

  File "/home/user/apps/anaconda3/lib/python3.6/xml/etree/ElementTree.py", line 1315, in XML
    return parser.close()

  File "<string>", line unknown
ParseError: no element found: line 1, column 0

As far as I can tell, I'm following the documentation correctly in my startup script, and my simple code roughly follows the GRASS wiki article on python.

What am I doing wrong? I've never used GRASS or GIS software before, so maybe I'm missing something obvious here. I'm using Grass 7.2 on a Debian Linux system with Python 3 (installed with Anaconda).


texts: Using OpenLayers and OpenStreetMaps to display statistical information with colors?

texts:
I'm planning to build a statistical tool using OpenLayers and OpenStreetMaps. This tool will have to display statistical data for French departments with colors. What I would like to do is the same as this :

http://sd1878.sivit.org/geoclip_stats_o3/index.php?profil=FR#v=map3;i=ocs.p_c1_2006;l=fr;z=178717,7038714,864706,751147

Does OpenLayers/OpenStreetMaps allow this ? Or should I look for other alternatives ?


texts:
I need to change the features of a MapLayer in real-time but I don't know how I can update/refresh/repaint the layer on the map without having to remove it and add it again (which creates flickering).

map.removeLayer(pathML);
SimpleFeatureCollection pointsFC = createFeatureCollection(points);
pathML = new MapLayer(pointsFC, pathStyle);
map.addLayer(pathML);

Do you know a better way to update the Layer rendering?


texts: Algorithm to fix topological errors in GIS software apps

texts:
Related question here.

I notice that ensuring topological correctness is essential for GIS applications, this is because the input from user or polygon boolean operation can have severe topological problems ( even though the polygons look correct) that would compromise on the quality of the subsequent operations.

Clean polygon is how Geo Wizards do to ensure the correctness of the topology.

Arcgis also has a command to clean up the slivers.

My question is not about how to use the existing software packages to ensure that the polygon input data is topologically correct; rather, my question is about how these software packages implement those cleaning procedures. In other words, what is the algorithm I can use to make sure that I can fix all the topological errors, given a set of polygonal inputs?


texts: Unable to remove join from two feature layers

texts:
Trying to remove the joins between two feature layers gets me this result after running each line individually through the interpreter:
MakeFeatureLayer("Y:/Desktop/Data/Levels_copy.shp","tempLevels")
<Result 'tempLevels'>
MakeFeatureLayer("Y:/Desktop/Data/Points_copy.shp","tempPoints" )
<Result 'tempPoints'>
AddJoin("tempPoints", "LEVEL_ID", "tempLevels", "LEVEL_ID")
<Result 'tempPoints'>
RemoveJoin("tempPoints", "tempLevels")
    Traceback (most recent call last):
      File "<stdin>", line 1, in <module>
      File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 6161, in RemoveJoin
        raise e
    ExecuteError: Failed to execute. Parameters are not valid.
    ERROR 000800: The value is not a member of  | Levels_copy.
    Failed to execute (RemoveJoin).

Am I missing something? Do the feature layers not exist in memory by the time I try to remove the joins?
Update
The joins still seem to exist as this codeblock still gives me fields from both shapefiles:
desc = arcpy.Describe("tempPoints").fields
for f in desc:
    print f.name


texts:
I have a PostgreSQL/PostGIS spatial relational database for my university project with objects, their coordinates and attributes (1:n, n:m relations).
I would like to map my spatial objects online and run queries, however, I am not ready yet to set up a web server and write server side scripts.

Are there any web mapping platforms that will allow me for now to visualize the data from the exported database? I suppose in this case I have to merge all of the related tables with attributes into one file.


texts: How to use Threads in PyQgis, mainly to keep UI active?

texts:
I'm trying to develop a plugin, which does heavy computations and as a result, user interface becomes inactive. I know that I should use threads to keep the UI alive (specifically, I have a progress bar that I want to update), but my experience with threads is very limited.

I was wondering if there are any good tutorials, or code snippets with documentation that cover this issue.
sample of my code so far (without threads):

    def initGui(self):
    # Create action that will start plugin configuration
    self.action = QAction(
        QIcon(":/plugins/lts/icon.png"),
        u"LTS calculator", self.iface.mainWindow())
    # connect the action to the run method
    self.action.triggered.connect(self.run)

    # Add toolbar button and menu item
    self.iface.addToolBarIcon(self.action)
    self.iface.addPluginToMenu(u"&LTS Calculator", self.action)

    # QtCore.QObject.connect(self.dlg.ui.find_cc, QtCore.SIGNAL("clicked()"), self.find_connected_components)

    QtCore.QObject.connect(self.dlg.ui.process_Button,QtCore.SIGNAL("clicked()"), self.process)
    # QtCore.QObject.connect(self.dlg.ui.layerCombo,QtCore.SIGNAL("currentIndexChanged(int)"), self.update_lts_field)
    QtCore.QObject.connect(self.dlg.ui.layerCombo,QtCore.SIGNAL("activated (int)"), self.update_lts_field)
    QtCore.QObject.connect(self.dlg.ui.find_cc_Button,QtCore.SIGNAL("clicked()"), self.find_connected_components)

    self.update_ui()
    self.layers = self.iface.legendInterface().layers()  # store the layer list
    # self.dlg.ui.layerCombo.clear()  # clear the combo
    for layer in self.layers:    # foreach layer in legend
        self.dlg.ui.layerCombo.addItem( layer.name() )    # add it t

   def process(self):
    """ Calculates Stress for the selected layer"""
    index = self.dlg.ui.layerCombo.currentIndex()
    if index < 0:
        # it may occur if there's no layer in the combo/legend
        pass
    else:
        layer = self.dlg.ui.layerCombo.itemData(index)
    # layer = QgsVectorLayer(self.fileName, "layer_name", "ogr")
    nFeat = layer.featureCount()
    layer.startEditing()
    ...
    street.compute_LTS()
        if street.LTS != None :
            i+=1
            j=ceil(i/(nFeat/100))
            self.dlg.ui.progress_bar.setValue(j)
            layer.updateFeature(feature)
# run method that performs all the real work
def run(self):
    # show the dialog
    self.dlg.show()
    self.dlg.ui.progress_bar.setValue(0)
    self.dlg.ui.progressBar.setValue(0)
    self.dlg.ui.layerCombo.clear()
    self.dlg.ui.lts_combo.clear()

    layers = QgsMapLayerRegistry.instance().mapLayers().values()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer and layer.geometryType() == QGis.Line:
            self.dlg.ui.layerCombo.addItem( layer.name(), layer )

    self.update_lts_field()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result == 1:
        # do something useful (delete the line containing pass and
        # substitute with your code)
        pass

Thanks


texts: Error getting spatial reference as string in arcpy

texts:
I am trying to get the string associated with a spatial reference to print into an xml using arcpy and xml.cElementTree.  Using the snippet:

...
spatial_ref = arcpy.GetParameterAsText(3)
...
Prj = et.SubElement(spaceref, "Prj")
space = arcpy.SpatialReference(spatial_ref)
Prj.text = arcpy.SpatialReference.exportToString(space)
...

Using essentially the same method, I am able to get the string to print with IDLE by pointing directly to a projection (.prj) file, e.g.:

  u"GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision"

but Arc returns this error when I try to run my script:

  Traceback (most recent call last):
    File "\DEFOE\Proc Storage\ORTHO\Arc Custom Tools\ImportUltracam.py", line 99, in
      space = arcpy.SpatialReference(spatial_ref)
    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\arcobjects\mixins.py", line 927, in init
      self._arc_object.createFromFile(item)
  RuntimeError: ERROR 999999: Error executing function.

Is there a way to input a spatial reference with Arc and return the string to be printed into the xml?

Update
Playing with IDLE, I tried with an EPSG code and was able to do exactly what I wanted:

>>> epsg = arcpy.SpatialReference(32618)
>>> print epsg
<geoprocessing spatial reference object object at 0x000000000C7D0710>
>>> arcpy.SpatialReference.exportToString(epsg)
u"PROJCS['WGS_1984_UTM_Zone_18N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-75.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]];-5120900 -9998100 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision"

So I tried with an EPSG code in the script:

spatial_ref = "32618"
Prj = et.SubElement(spaceref, "Prj")
space = arcpy.SpatialReference(spatial_ref)
Prj.text = arcpy.SpatialReference.exportToString(space)

And I get the same error.


texts: Analyzing large rasters from an image service

texts:
I'm working with NAIP imagery in ArcGIS 10.4.1 that I'm accessing through the USDA server (https://gis.apfo.usda.gov/arcgis/services).

I want to perform focal statistics on bands 2 and 4 (green and IR), which is easy when the file is hosted locally because I can just double click on the raster and select any of the individual bands.

With this NAIP file however, navigating to the hosted folder in the focal statistics window shows an empty folder. Clipping the NAIP down and saving as a local file is time-consuming and seems to warp the display of the image.

Is it possible to extract bands from non-local rasters?


texts: Generating consistently-dimensioned polygons in mm units?

texts:
I have a function which creates Solar Photovolatic panels represented as
polygons. Essentially, it creates a rectangular grid where the user can specify the
following parameters:

Length
Width
Horizontal distance
Vertical distance

The code is based on the plugin FeatureGridCreator but focusing only
on the polygon aspect. It works well for the most part, especially when creating
polygons with large dimensions (e.g. 10m length and width; 10m horizontal and vertical distances).

But I noticed a couple of issues:

When specifying polygons for dimensions less than 2m for both length and width, no polygons were created.
When specifying polygons with differing dimensions (e.g. 5m length and 7m width), the dimensions were not the same when measured with the Measure Line tool. For these dimensions, the length and width were shown to be 4m and 6m respectively.

The CRS used for both the projection and layer is EPSG:27700 although I wouldn't have thought this would be a problem.

So does anyone have any idea what could be causing these issues? I am also open to suggestions as to how the code could be improved or even replaced with a better alternative.

Here is the code which can be reproduced in the Python Console, a polygon layer must be selected with a relevant CRS before running the function:

from PyQt4.QtCore import QVariant
from math import ceil

def generate_pv_panels(length, width, distance_x, distance_y):
    # Define layer properties
    layer = iface.activeLayer()
    crs = layer.crs()
    memory_lyr = QgsVectorLayer("Polygon?crs=epsg:" + unicode(crs.postgisSrid()) + "&index=yes", "PV panels for " + str(layer.name()), "memory")
    QgsMapLayerRegistry.instance().addMapLayer(memory_lyr)
    memory_lyr.startEditing()
    provider = memory_lyr.dataProvider()
    provider.addAttributes([QgsField("ID", QVariant.Int)])
    fid = 0
    start_x = 0
    start_y = 0
    # Ensure polygons are not created 'within each other'
    if distance_x < (length / 1000):
        distance_x = (length / 1000)
    if distance_y < (width / 1000):
        distance_y = (width / 1000)
    fts = []
    for f in layer.getFeatures():
        fid += 1
        bbox = f.geometry().boundingBox()
        start_x = bbox.xMinimum() + float(distance_x / 2)
        start_y = bbox.yMinimum() + float(distance_y / 2)
        for row in range(0, int(ceil(bbox.height() / distance_y))):
            for column in range(0, int(ceil(bbox.width() / distance_x))):
                fet = QgsFeature()
                geom_type = pv_panel_size(length, width, start_x, start_y)
                if f.geometry().contains(geom_type):
                    fet.setGeometry(geom_type)
                    fet.setAttributes([fid])
                    fts.append(fet)
                start_x += distance_x + (length / 1000)
            start_x = bbox.xMinimum() + float(distance_x / 2)
            start_y += distance_y + (width / 1000)
    provider.addFeatures(fts)
    memory_lyr.updateFields()
    memory_lyr.commitChanges()

def pv_panel_size(length, width, x, y):
    # Length & width measured in mm; x & y measured in m
    l = length / 2000
    w = width / 2000
    return QgsGeometry.fromRect(QgsRectangle(x - l, y - w, x + l, y + w))

generate_pv_panels(10000, 10000, 100, 100)


texts: Using ModelBuilder and Raster Calculator with Output Location?

texts:
I'm using ArcMap 10.3.1 ModelBuilder to create a tool that calculates the surface temperature using with base a landsat 8 bands. The tool is working fine, my parameters are working well and inputs too.

My only problem is that I want to the system be most automated as I can, and I want the user to be capable to set the output location of the files and its name. Using another tool's at ModelBuilder I checked that is easy to add it as model parameter, but at raster calculator these settings are not available.

Does somebody know how to make available to the user to set the output file location, and name using the ModelBuilder and raster calculator tool?


texts: Auto generating waypoints to a GPX file along track (hiking path - not as crow flies)

texts:
I have used Google Maps to build a list of waypoints for a trail I want to hike but the waypoints are simply my start, end and maybe a couple of points of interest for my day and any campgrounds along the route.

I now have all of this data in a GPX file for importing into my GPS, but it keeps reporting distance as the crow flies to the next waypoint instead of distance along the trail. I have since learned that the fix is to add a metric boatload of waypoints to my tracks.

Is there any way to auto-generate waypoints along the track line between two existing waypoints that follows the trail line instead of line-of-site?

I keep seeing people pointing to GPSBabel (but that only does line-of-site) or suggestions to manually add points.


texts: GeoPandas: Get distance from geom in tableA to a geom in tableB that matches attribute from tableA

texts:
UPDATE WITH SOLUTION (at bottom)

Using GeoPandas I am trying to identify the distance between an election candidate's home address (Point) and the electoral district they are running in (Polygon). Often they don't live in their district and so I'm trying to calculate how far each candidate lives from their riding as a crow flies. (I am unsure at this point if GeoPanda's distance method calculates to the nearest vertex/node or the centroid when using a polygon, but at this point its not important)

I have a DISTRICT GeoDataFrame which contains the electoral district boundaries:

id  ed_name     geom
0   MapleRidge  POLYGON(599240.6488817427....)
1   St.Johns    POLYGON(589240.6488427823....)
2   Southgate   POLYGON(563405.6488424563....)
3   etc....     etc...

And my second GeoDataFrame -- CANDIDATES -- contains details about the candidates, including the POINT geom of their home address and the district they are running in.

id  candidate  ed_running_in   geom
0   John       MapleRidge      POINT (523640.6482456427....)
1   Steve      Southgate       POINT(659240.6588817427....)
2   Shelly     St.Johns        POINT(879240.6488817427....)
3   Irene      MapleRidge      POINT(129240.6288817427....)
4   Patrick    MapleRidge      POINT(659240.1688817427....)
5   Ian        Southgate       POINT(929240.9888817427....)
6   etc...     etc...          etc...

(Note that these are fictional tables for illustrative purposes, so the geoms are bogus.)

GeoPandas has the following method: GeoSeries.distance(other).

I want to add a columns to the CANDIDATE table called "Distance" that shows how far the candidate lives from the riding they are running in. The Distance method in geopandas calculates element wise which is not useful here. I need the distance method to choose the correct district to measure against candidate home based on the "ed_running_in" value from the CANDIDATE table.

In PostGIS this is easy by using something to the effect of:

WHERE candidate.ed_running_in = district.ed_name

Below is some invalid code (loc doesn't allow you to reference another df)  that illustrates what I'm mechanically trying to accomplish after setting the "ed_name" field as the index for the DISTRICT table:

candidates['distance'] = candidates['geom'].distance(district.loc[candidates['ed_running_in']])

UPDATE AND SOLUTION:
The key is to build a temporary bridge table, which merges over the geometry from the DISTRICT GeoDataFrame (i.e. the boundary file) to a cloned CANDIDATES table.
Because the distance method is comparing two geometries pairwise you must create a bridge table of equal length and in matching row order.

So you have to clone the CANDIDATES table (say: CANDIDATES_TEMP) and then merge the geometry from districts over to it based on a match with "ed_running_in" == "ed_name".
But the key thing (which took me a while to figure out), is that you must first drop the geometry column (the home addresses point) from the from cloned table, because GeoPandas will not accept two geometry values in a single element -- if you do not do this, it will convert your geometries to a string and therefore also converts the GeoDataFrame to a simple DataFrame. That is really the core Eureka moment here.

From here, you simply have to use the GeoSeries.distance(other) method to get the distance between the POINTS geometry from the CANDIDATE table to the POLYGON geometry from the cloned CANDIDATES_TEMP table.

PS. And in case you're wondering about the behaviour of the distance method when measuring a POINT against a POLYGON: the calculation is based on the nearest vertex, and not the centroid.


texts: Transformation from CRS EPSG:3857 to ESRI:54009

texts:
I am trying to convert from EPSG:3857 to ESRI:54009 and facing some issues (pyproj version -2.6.1).
ERROR 1: PROJ: proj_create_from_database: crs not found
ERROR 1: PROJ: proj_create: unrecognized format / unknown name
ERROR 6: Cannot find coordinate operations from EPSG:3857' to `'``
My code as following:
def project_geometry(geom , in_crs = 3857 ,out_crs = 54009):
    inSpatialRef = osr.SpatialReference()
    inSpatialRef.ImportFromEPSG(in_crs)
    outSpatialRef = osr.SpatialReference()
    outSpatialRef.ImportFromEPSG(out_crs)
    coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
    geom =ogr.CreateGeometryFromWkb(geom.to_wkb())
    geom.Transform(coordTrans)
    return loads(geom.ExportToWkb())

Is there a direct way to do this transformation?
or do I need to use transformation from EPSG:3857 to other CRS and than from that to ESRI:54009?
Also I didn't find transformation files in this pyproj release.


texts:
I am trying to get my labels transparent. I have converted them into annotation and unable to make the annotations transparent. I have kept annotation in the map. My annotation layer does not appear on table of contents!
Ideas of getting the lables transparent or annotations transparent will be appreciated.


texts: QGIS QgsAttributeTableModel only loads when called without a function in Python

texts:
I tried to load QgsAttributeTableModel into QgsAttributeTableView or a plain Qt QTableView but it could run from inside a class method or function but it runs without any function or method. Below is the code that doesn't work.

This one works...

from PyQt4.QtGui import QTableView, QDialog, QVBoxLayout
from qgis.core import QgsVectorLayerCache
from qgis.gui import (QgsAttributeTableModel,
                      QgsAttributeTableView,
                      QgsAttributeTableFilterModel)
import qgis.utils
layer = iface.activeLayer()
if not layer is None:
    canvas = iface.mapCanvas()
    vector_layer_cache = QgsVectorLayerCache(layer, 10000)
    attribute_table_model = QgsAttributeTableModel(vector_layer_cache)
    attribute_table_model.loadLayer()

    attribute_table_filter_model = QgsAttributeTableFilterModel(
        canvas,
        attribute_table_model
    )
    attribute_table_view = QgsAttributeTableView()
    attribute_table_view.setModel(attribute_table_filter_model)
    dialog = QDialog(iface.mainWindow())
    attribute_table_view.setParent(dialog)
    layout = QVBoxLayout(dialog)
    layout.addWidget(attribute_table_view)
    dialog.resize(800,600)
    dialog.show()

But this one doesn't

from PyQt4.QtGui import QTableView, QDialog, QVBoxLayout
from qgis.core import QgsVectorLayerCache
from qgis.gui import (QgsAttributeTableModel,
                      QgsAttributeTableView,
                      QgsAttributeTableFilterModel)
import qgis.utils

def open_attr_table(layer):
    if not layer is None:
        canvas = iface.mapCanvas()
        vector_layer_cache = QgsVectorLayerCache(layer, 10000)
        attribute_table_model = QgsAttributeTableModel(vector_layer_cache)
        attribute_table_model.loadLayer()

        attribute_table_filter_model = QgsAttributeTableFilterModel(
            canvas,
            attribute_table_model
        )
        attribute_table_view = QgsAttributeTableView()

        attribute_table_view.setModel(attribute_table_filter_model)
        dialog = QDialog(iface.mainWindow())
        attribute_table_view.setParent(dialog)
        layout = QVBoxLayout(dialog)
        layout.addWidget(attribute_table_view)
        dialog.resize(800,600)
        dialog.show()
layer = iface.activeLayer()
open_attr_table(layer)

Both codes are identical. The only difference is, one runs inside a function, the other doesn't. I have tested them in QGIS 2.8 and 2.14.

I am wondering if QGIS has disabled the use of QgsAttributeTableFilterModel and QgsAttributeTableModel in QGIS plugins. Is there anything I should do to fix this issue?

Similar issues with no solution

show-attribute-table-in-the-form-table-of-contents

how-to-show-attribute-table-in-the-form


texts: How to make a shapefile (which is not grid aligned) to grid aligned in QGIS?

texts: I have shapefile which is not grid aligned. I want it to be grid aligned. Doing it manually will consume lot of time. Is there any better way to make the shapefile to grid aligned?


texts: Is there a recomended method for symbolizing inland lakes, streams and rivers in a topographic map?

texts:
I produce a lot of maps with inland lakes, rivers and streams overlaid over a dem surface symbolized with a color graduation (green for lows and brown for higher elevations) I also use a semi-transparent shaded relief to enhance the terrain. I have a standard blue color I often used for both lakes and streams but I noticed that sometimes the hydrology does not stand out enough. I experimented with outline for the lakes and with dark and light blues but nothing seems to work for every map.

Can someone recommend a good color/method for symbolizing lakes/streams/rivers?

Is there a generally accepted standard theme for symbolizing lakes/streams/rivers? (Ontario, Canada)


texts: What do "count", "mean", and "sum" in QgsZonalStatistics mean?

texts:
I am working with Qgis and making some analytics of lightdata.
I've used these commands:

rasterfile = qgis.utils.iface.mapCanvas().currentLayer().source()
vectorlayer = qgis.utils.iface.mapCanvas().currentLayer()
import qgis.analysis
zonalstats = qgis.analysis.QgsZonalStatistics(vectorlayer,rasterfile)
zonalstats.calculateStatistics(None)

When I open the table I get mean, count, sum.
What does it mean in this term? count - what is it counting? sum - what is it summarizing? mean - how does it calculate the mean?


texts: Making cumulative sum of daily precipitation raster in ArcGIS Desktop?

texts:
I have daily precipitation rasters from 1981 to 2016 and i need to get the total sum of each pixel in the raster and the mean of them.

How can i iterate through them to sum all of them (more than 13000 rasters)?

I've tried Cell statistics with the iterate rasters loop in arcmap but it does not work.


texts: How to put geotagged .jpg thumbnails on a map in QGIS?

texts:
I have a collection of photos (already geotagged in EXIF) and I'd like to put thumbnails of these pictures on a map (at the location where they were taken).
Just looking for a static visualisation of a map with picture thumbnails - no need for any interactivity (the thumbnails need not to be clickable or mouse-overable), but I do need a way to control the size of the thumbnails (but it can be a global setting that applies to all the pictures). There would be only a handful of thumbnails on a map.
I am new to QGIS, and hoped that I could accomplish this by creating raster layers from .jpg images - and although the pictures are read into QGIS and the coordinates recognised (in Properties), I don't see the picure in the right place and I cannot resize them.
I tried looking into eVis, but that also doesn't seem to do what I want.
Oddly enough, also in Google I cannot find any examples of QGIS-made maps that have these .jpg thumbnails on them. With all the dazzling options in QGIS, surely there must be a simple way to achieve this?
Edit: While still researching this, I came across the digiKam application (digikam.org) that does exactly what I would like to do in QGIS. It puts photos on a map layer, and lets you change the size of thumbnails (see screenshot) Unfortunately, the thing it doesn't do is read GPX files and other GIS functionality...


texts: APDM - error exporting database from Oracle

texts:
I have a fairly large database on a slow network and I am using a third party tool to update data which calculates other fields and so on.

I have a guid string with an empty box (weird glyph - seen in shapefile tables before) in the middle of the string for just a few records, seemingly random.

{D62D[][]1F1-CA71-4490-B940-046E61989DA6}
{A16659DF-7131-4[][]E-8DFF-DBF71257DC08}

The empty boxes are not brackets. I just did that for the example - they are empty boxes.

I would like to use python or some other script/sql to run the entire database 300,000 records to see if this is present elsewhere?
What are possible causes?
How would this be related to exporting the database?


texts: Displaying map labels over coordinate grid in QGIS?

texts:
My grid overlaps labels placed on a map.

Is there a way of forcing the labels to be displayed in front/on top and the grid placed behind the text?

Despite looking I cannot find a solution.

Using QGIS 2.18.15


texts: Error Using OpenLayers 3 Web mapping applications

texts:
As being new to web site development, I have been able to display map on my website, but I want that I should only zoom to some specific places of India on loading of web pages.

Can somebody please suggest on how to go about it.

<!doctype html>

<meta charset="utf-8">
<title>Web Urban Atlas</title>
<link type= text/css rel = stylesheet href="urban_stylesheet.css" />
<link rel="stylesheet" href="http://code.jquery.com/ui/1.11.2/themes/smoothness/jquery-ui.css">

Urban Atlas

Bangalore

Bangalore is city of gardens

delhi

Content for Delhi

mumbai

Content for Mumbai


texts: Multi-writer (ArcGIS FME) for export to DGN

texts:
I would like to ask, if there is any way for export from one Reader (line geometry, format ArcSDE) to more writers (all in the same format DGN). Specific line with identical ID then would be displayed in Microstation for example twice, but with different symbology - lower line will be wider with different color and upper line will be thinner. I would like to induce effect of line´s border.

It would work together a connection of another writers, which would export to another formats (database) of output datasets except DGN?


texts: QuickMap Services plugin for QGIS corrupted

texts:
I probably need to use QuickMapServices plugin but when I try to install it I receive thi error message:

I use QGIS 2.6.1 on windows 10 64bit

I am not the admin of the PC so I am not able to update to newset version of QGIS. The plugin page shows that the minimun version is 2.4.0


texts: How to create a real size buffer around a point in qgis

texts:
I have a postgis pointlayer in EPSG:3857, added to a project that is set up with OTF ON to project everything (other layers) to EPSG:3857.
In postgis, I defined a view to dynamically create a buffer around the points, with a radius defined by an attribute:

CREATE VIEW point_buffer AS
     SELECT gid,st_buffer(the_geom,buf_radius) AS GEOM
FROM points_layer

It creates the buffer correctly, but when I measure the radius of the buffer in qgis, it shows the buffer without taking the ellipsoid transformation into account.
(when i measure the buffer radius with OTF OFF, it shows the correct radius)

Is there a way to have a buffer around a point be scaled to the ellipsoid?


texts: How to Create Spatialite tables and layers in QGIS using Python

texts:
I'm trying to build a data model for use in QGIS. I want to use Spatialite (SL). I've found a few things on the internet to create an in memory layer. I don't know how to specify the SL database where I want it to go. I also stumbled across this article, and I wasn't sure if this was still valid (over 2 years old).
How to create a SpatiaLite Layer in QGIS 3 with Python?
Below is the code that I found to create the layer, but I can't find how to create it in SL.
uri = "point?crs=epsg:4326&field=id:integer"
scratchLayer = QgsVectorLayer(uri, "Scratch point layer",  "memory")
TYIA


texts: how to draw curved line in a qgsmapcanvas

texts:
i am working on a project with qt4.8.6, api c++ of qgis 2.8.2
we need to draw curved line in a vectorlayer by given 2 points and an offset by exemple

is it possible to do that with the c++ api of qgis ?


texts: Google Earth Engine - Name Layers with Image Image ID

texts:
I have recently started to explore google earth engine - and would like to know if there is a way to name image layers (within a collection) with the individual image ID? I have come this far:

var s2 = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(roi)
  .filterDate('2019-03-01','2019-03-15')
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
  .limit(5);

  var rgbVis = {
  min: 0.0,
  max: 3000,
  bands: ['B4', 'B3', 'B2'],
};

function addImage(image) { // display each image in collection
  var id = image.id;
  var iMage = ee.Image(image.id);
  Map.addLayer(iMage, rgbVis, 'RGB');
}

s2.evaluate(function(s2) {  // use map on client-side
  s2.features.map(addImage);
});


texts:
I have problems when I want to log the administrative user through GeoExplorer when geoserver is running on a server with public ip but redirection to a private ip.

GeoExplorer takes too long to load and when it finally starts the login button does not appear in the interface. (in firedebug can be seen that there is no of cookies data).

If I access the GeoExplorer from the internal network, the login button load properly, GeoExplorer fast loading and allows me to make the user login.

The only open through public ip port is 80. Instead for access through the internal network are all open ports.

I am using: geoserver 2.4.0 , Ubuntu 12.04 LTS


texts: Printing centered map from a QGIS project for each point in shapefile?

texts:
I need to produce on the order of 100 maps that are centered on each point of interest in a shapefile. I would like to prepare all the layers in a master QGIS project, and set up the composition for one point (so that printing 100 maps could be done manually, if need be).

I'd have something like the following layers:

basemap
points
point buffer

And I would like to then automate printing to svg something like:

For each point in a shapefile
Center the map canvas on that point
Filter a buffer layer to only have that point
Print an svg with composition from the map composer

I'm reasonably certain I know how to do 1 & 2, but haven't found details on 3 & 4 on this site.


texts: ArcGIS online Hosted Feature Service will not Open ArcGIS Desktop if only shared with a group

texts: I have an ArcGIS Online Hosted Feature Service shared with a group. When I try to open that feature service in ArcGIS desktop the feature service will not open. When I change the sharing capabilities to public the feature service will open in ArcGIS desktop. The group was created after July 2015 and I am the group administrator. Also I am signed into my ArcGIS online account when I have ArcGIS desktop open. I have also added my organizations URL in ArcGIS Administrator, and I am able to publish Feature and Map services.


texts: Clip Shapefiles via Iteration in QGIS Python Script

texts:
I have two polygon shapefiles with several objects. Each object in both shapefiles has a given ID number. Now I want to clip the object with ID=1 in Shapefile A with the object with ID=1 in Shapefile B.

However each shapefile has about 5000 objects, so I would like to create a Python script via the QGIS graphical modeler. I know how to use the graphical modeler to clip objects and stuff, but I cannot find any iteration option within the modeler. Is there an option to do so or how can I implement an iteration directly into the script to get the wanted clip results?

In the end I would like to have one new shapefile with all the clips.


texts: Check if feature is selected (Openlayers)

texts: I wonder if there is a way to check whether a pointfeature (OpenLayers.Feature.Vector) is selected. I read the API, but didn't find any solution. Any idea?


texts: geopandas.to_crs(epsg=3395) does not work in jupyter notebook. it seems some mismatch with pyproj & geopandas versions?

texts:
After importing geopandas df , with epsg =4326, I tried to re project it to epsg =3395 using
gdf.to_crs(epsg =3395)  in jupyter python note book.
I got error "RuntimeError: b'no arguments in initialization list' "
I have ,
geopandas 0.6.3
pyproj 1.9.6
but I could not update to latest version using "conda -forge channel"
I also have proj4 5.2.0 in base channel  , should I remove it ?


texts: CartoDB Torque animation running longer than 'duration' value in CartoDB editor

texts:
I'm creating a time series Torque animation with CartoDB, using data that spans a one week period. I've noticed that the animation is running longer than the value that 'duration' is set to in the CartoDB editor. This is true no matter what value I use for the duration. For example, if I set duration to 30 seconds, the animation actually takes about 35 seconds to complete one full cycle. If I set it to 15 seconds, it seems to take about 17.5 seconds.

The map and data can be viewed here.

The relevant CartoCSS code is

Map {
-torque-frame-count:1024;
-torque-animation-duration:30;
-torque-time-attribute:"timestamp";
-torque-aggregation-function:"count(cartodb_id)";
-torque-resolution:2;
-torque-data-aggregation:linear;
}

Does anyone have experience with this issue or thoughts on what might be causing it?


texts: Where can I find a simple TOWNSHIP/RANGE, SECTION, QUARTER SECTION shapefile for Adams County Colorado?


texts: Continuously Updating Operations Dashboard Without Breaking Widgets?

texts:
I am a new user of ESRI's Operations Dashboard and I am experiencing a problem. I am working on a continuously evolving project, so the data in my attribute table changes each day and gets updated every evening. What I am looking for is a way to bring in the new update attribute table data every evening and have my widgets automatically connect to it and work correctly. As of now when I update my attribute table data and republish it online it breaks all of the widget and I have to recreate them all.

If I had a "Data" layer that my widgets were pulling from and I overwrote it with a new layer named "Data" would that fix my problem?

Just trying to save myself from recreating my widgets every day.


texts: Calculating st_simplify tolerance for OpenLayers

texts:
I want to display a vector layer in OpenLayers. Users can scale the map before they load the vector layer and I'd like to simplify the geometry with st_simplify before sending it to the client.

Does anyone have a simple formula I can use to calculate an appropriate st_simplify tolerance based on the map's properties (e.g. zoom, scale or extent?)?

https://postgis.net/docs/ST_Simplify.html
http://dev.openlayers.org/docs/files/OpenLayers/Map-js.html

(I don't think this is relevant but have a PHP middle layer and am passing the layer as KML)

I've searched around and maybe I'm looking the wrong way because I can't believe I'm the first one with this question!


texts:
I've a MODIS04 L2, and I want to extract band
8:Image_Optical_Depth_Land_And_Ocean. I use Quantum GIS to open this band,
because Quantum GIS see this is not coordinates so this allow me to select
projection as below: +proj=utm +zone=50 +ellps=WGS84 +datum=WGS84 +units=m
+no_defs +towgs84=0,0,0.

All I know is that I want to translate to WGS84. This open this band seems
good but when I open this band with ENVI this is better (differents output
images). But I see the Modis Conversion Toolkit use same warp parameters as
Quantum Gis expect it has a parameters ("Numbers of Warp Points").

I think that I miss something with QGIS and GDAL too.

This is Envi Conversion Toolkit

This is open file HDF with QGIS


texts: Merging two image collection into one with two bands (LST and NDVI )

texts:
I have created two image collections (monthly average) from different datasets, one is LST (MODIS MOD11A2 ) and the another one is NDVI (MODIS MOD13Q1). The two datasets have exactly the same date.

How to merge them together into one image collection, which each image has two bands (LST and NDVI) ?


texts: How to visualize multiband imagery using rasterio?

texts: How do I plot a satellite image using rasterio? I'm able to plot the individual bands but not the entire image. Is there a way to combine the separate R, G, B numpy arrays to create a standard RGB image?


texts: Clustering polygon's created with Leaflet.Draw

texts:
I am experimenting with Leaflet at the moment and am trying to cluster polygon's that I've created with Leaflet.Draw plugin's tools. I know that clustering polygons is possible (Is it possible to cluster polygons in Leaflet?) and have successfully clustered polygons that I've created manually by creating an extended ClusterablePolygon class. I can also cluster markers that are created by Leaflet.Draw. This brings a twofold question:

Is it possible to alter the draw tools in such way that it automatically creates a custom ClusterableGeometry instead of a regular geometry?
Is it possible to instantly after creation convert a regular geometry into ClusterableGeometry?


texts:
As part of a larger program (basically a remake and extension of DistanceCalc) I want to programmatically create a line between the two closest points of two regions in MapInfo.

The ObjectDistance function effectively measures this line and returns its length but there is no way to store the line itself.

The "Nearest" statement almost does what I want but it does too much! - it measures from one site (or each item in a table) to each item in another table. Also it puts the results in a table whereas I want to manipulate them before storing them.

Is there a simple way to do this?


texts: Convert latitude longitude coordinate to x y coordinate?

texts:
I'm using a postgis database and java to develop my application.
I need to transform my geographic data currently in latitude/longitude format into a x/y format.
I found solutions based on UTM coordinates, however points in different zones are a problem.
I'm looking for a java library or a algorithm that does exactly this:
http://www.whoi.edu/marine/ndsf/cgi-bin/NDSFutility.cgi?form=0&from=LatLon&to=XY

Given two lat/lon coordinates, one of them for origin, the output is x and y in meters.
Thank you for your time


texts: Field calculator to copy numeric values from other fields

texts:
FieldA contains text: both numeric(123) and alpha-numeric(123a, 2-3, 3/2, '). FieldB is longstring. I want to use field calculator and codeblock  on FieldB to bring from Field A only the numeric values.

while True:

try:
    val = int(x)
    print (val)
    break
except ValueError:
    print ("non-numeric")
    break

I tried something similar to the code seen here, but the process fails at some point. I just need something to see if fieldA .isdigit and if true copy to fieldB(if false pass a string "non-numeric"). Or anything similar.

I am using ArcGIS 10.2 for Desktop.


texts: Shapefile calculation error: 'list' object cannot be coerced to type 'integer'

texts:
I'm trying to identify which of two regions within a shapefile specific spatial points lie. To this I used this code.

This worked well but I ended with a lot of NAs because my 1/3 of my points are slightly outside the shapefile boundaries.

I ended up creating a new shape file with a buffer around the two regions. However, for some reason, despite the shape files have the same class columns as the old shapefile etc, it does not run and I get this error

map = read_sf("OSPAR23_buf1.shp")
pnts <- subset(SENT, select=c(Longitude, Latitude, ID))

pnts_sf <- st_as_sf(pnts, coords = c('Longitude', 'Latitude'), crs = st_crs(map))

pnts <- pnts_sf %>% mutate(
  intersection = as.integer(st_intersects(geometry, map)),
  OSPAR = if_else(is.na(intersection), '', map$Region[intersection]))

Error: 'list' object cannot be coerced to type 'integer'

I then tried to create a buffer around the pnts file using this code taken from this post

pnts_sf_buf = st_buffer(pnts_sf, 0.25)

And run the previous code with the maps that had worked previously. But I get the same error.

There does not seem to be any obvious difference between the shape files that worked and those that don't, and I really can't figure out why this isn't working.

A subset of my data is below. A zipfile of the old shapefile that worked is here. The new shapefile that gives the error is here

structure(list(Longitude = c(-3.55322138, -4.694794872, -6.116016097,
-4.520111046, -0.910337034, -3.408259063), Latitude = c(50.43247257,
50.34775167, 57.13215105, 51.73420664, 50.78851095, 54.27026104
), ID = 1:6), row.names = c(1L, 2L, 3L, 6L, 7L, 14L), class = "data.frame")


texts: ArcGIS : You do not have permissions to use this resource. (Error Code: 401)

texts:
I am trying to query Demographic data using ArcGIS Python API and the following is my code to do the same.

from arcgis.gis import GIS
from arcgis.geoenrichment import Country, enrich, BufferStudyArea

gis = GIS("https://www.arcgis.com",myusername, mypass)

prt = Country.get('PRT')
df_prt = prt.data_collections

enrich(study_areas=[{'address':{'City':'Mafra'}}], data_collections=['PopulationTotals'])

And, on running the last line, I get the following error :
Exception                                 Traceback (most recent call last)
<ipython-input-27-f4e6881f3364> in <module>()
----> 1 enrich(study_areas=[{'address':{'City':'Mafra'}}], data_collections=['PopulationTotals'])

4 frames
/usr/local/lib/python3.7/dist-packages/arcgis/gis/_impl/_con/_connection.py in _handle_json_error(self, error, errorcode)
    534
    535         errormessage = errormessage + "\n(Error Code: " + str(errorcode) +")"
--> 536         raise Exception(errormessage)
    537     #----------------------------------------------------------------------
    538     def post(self,

Exception: You do not have permissions to use this resource.
(Error Code: 401)

I am using ArcGIS for the first time, and I'm not sure what is causing the problem. I only followed the documentation.
What is the problem and how can I get rid of this error?


texts: Apply function from rgeos package to data sets

texts:
I am using gDistance function from rgeos package to calculate the nearest distance from points to a polygon.
Originally I used for loop to calculate the distance for each row in my points data set, however, the speed is extremely slow, especially when I have a large data set. So I tried to use 'apply'.
My codes look like:

Invasive_dist <- apply(sp_invasive_pixels, 1, function(x) gDistance(x, sp_tumor))

whereas sp_invasive_pixels is my SpatialPoint object, and sp_tumor is the SpatialPolygons object.
My code does not work, the error message is:

  Error in apply(sp_invasive_pixels, 1, function(x) gDistance(x, sp_tumor)) :
    dim(X) must have a positive length

What do I need to change?

test program:

test_poly <- data.frame(matrix(nrow = 4, ncol = 2))
test_pts <- data.frame(matrix(nrow = 4, ncol = 2))

test_poly$X1 <- c(1, 1, 2, 2)
test_poly$X2 <- c(1, 2, 1, 2)

test_pts$X1 <- c(4, 5, 6, 7)
test_pts$X2 <- c(5, 6, 7, 8)

sp_pts <- SpatialPoints(test_pts)

sp_poly <- SpatialPolygons(list(Polygons(list(Polygon(test_poly)),1)))


texts: Plot points with coordinates from txt-file

texts:
I'm using QGIS 3.12.

I have a list of coordinates in txt file, as follows

"longitude" "latitude"
48.819056, 2.127009
48.782213, 2.143340
48.801580, 2.355238
48.904471, 2.296038

How to plot this list of points on a layer in my project?


texts: Can Compound Topographic Index (CTI) / Topographic Wetness Index (TWI) produce negative values?

texts:
I am calculating Topographic Wetness Index for a study area using SRTM DEM (90 and 30m).

I've seen several versions of TWI (CTI).
I have been using: TWI = ln(As +1 / tanB)
      where As = catchment area, B = slope (degrees/radians??)

When I run this function, the output contains negative values (using slope in degrees and radians).
Are negative values acceptable? Or does this indicate an error in the analysis??

I understand the output should be scaled to 1-10, but negative values seems counter intuitive. Could someone explain what negative values could mean (due to flat areas, no data values?)

On a side not, the calculations are being run on on a very large area, which overaps several watersheds. Is it appropriate to run TWI over such a large extent?


texts: Why Attribute name length in QGIS is limited to 10 characters?

texts:
In QGIS software, a vector layer has "attribute" but the name length is limited to only 10 Characters.

Is there any method to name an "attribute" to a longer name? (e.g. attribute name "building name" becomes "building n")

My layer define dialog window is as below


texts: Google Map layer loads tiles even if this layer is disabled in OpenLayers

texts: Sorry for cross posting I've asked this question in OL mailing list, but not get response yet. It is very important for me to resolve this issue. I've made a little example that reflect my problem. It's not clear for me why GM layer loads tiles even if this layer is disabled (it can be observed against firebug in Firefox) and how I can to prevent this?


texts: how to use a database's table data in QGIS

texts: I have already loaded a MySQL database as a layer, one of the tables contains data that are coordinates of some places, how to plot these data on a map(for example google maps)?


texts: Export a style from Mapbox Studio to be used in Mabox GL-JS

texts:
I designed my map with Mapbox Studio.

However I would like to use my design in Mapbox gl-js. The styles are not CartoCSS but a JSON.
Is there a way to export my design to this JSON format?

They seem very close, so I guess I could write a converter, but I would require some time, and I didn’t find all the features in the JSON format (for instance selecting with a regex).


texts: Converting dimension feature into image using ArcObjects?

texts:
I have a task to convert dimension feature into image with ArcObjects.

I need to get real image that I see in ArcMap (with right size and spatial reference, of course) instead of turning feature into raster by values of some field.

Is it possible to perform feature-to-image conversion with ArcObjects?

I need solution applicable to just dimension feature class, not to layer on the map.

It seems that the IDimensionGraphic interface should be used. But I don't understand what arguments should be passed to its Draw method.


texts:
Is it possible to load GeoJSON directly into qgis from github?

The GeoJSON is accessible here https://raw.githubusercontent.com/garynobles/data/master/sites.geojson

It displays fine if downloaded as a local file but I would like to be able to access it directly without having to download. I have googled around, but everything is about exporting to GeoJSON for github (i.e. keeping crs wgs84).


texts:
I want to convert a GeoPandas GeoDataframe to raster (.tif).
I have multiple attributes in my GeoDataframe and want to create a raster using  a single field from the attribute table. The underlying vector format are 16x16m square polygons, which should be represented later by a single pixel.
Up to now I haven't found any straight-forward solutions or useful tutorials, so I am looking for some helpful answers or hints of where to find practical information.
I'm using Windows 10 and Python 3.8.


texts: ArcGIS Server - Print Map aynschronously using Map_Only template

texts:
I'm trying to print a Map asynchronously using the ArcGIS printing tools. I'm successfully able to print the map using any template other than "MAP_ONLY". When i use "MAP_ONLY" template i get this error :
ERROR 001305: Missing or invalid 'exportOptions' property in WebMap; cannot derive a size for the output image. Failed to execute (Export Web Map). Failed to execute (Export Web Map Task). I'm not sure if there is a way to include the export options in the web map Json object. Here is the code:

 var initPrintParameters = function() {
        var temp = new PrintTemplate();
        temp.exportOptions = {
          width: dpi * 10, //3000, //960, //840,
          height: (20 * dpi) / 3, //640, //560,
          dpi: dpi
        };

        temp.format = "PNG32";
        temp.layout = "MAP_ONLY";
        temp.preserveScale = true;
        temp.showAttribution = true;

        printParameters = new PrintParameters();
        printParameters.format = "png32";
        printParameters.template = temp;
        printParameters.map = map;
        printParameters.outSpatialReference = map.spatialReference;
        printParameters.async=true;
    };

  var exportMapAtCurrentExtent = function() {
    printTask = PrintTask(printTaskUrl);
    initPrintParameters();
     var jsonMap = JSON.toJson(printTask._getPrintDefinition(map,printParameters));
    var printLayout={
      Web_Map_as_JSON: jsonMap,
      Format: "PNG32",
      Layout_Template: "A3 Landscape"
    };
      gp.submitJob(printLayout,gpJobComplete, gpJobStatus,gpJobFailed);
  };

  var gpJobComplete=function(jobInfo){
      console.log("Job Complete");
      console.log(jobInfo);
      gp.getResultData(jobInfo.jobId,"Output_File",getResultUrl);
  };

  var getResultUrl=function(result){
      console.log(result.value.url);
  };

  var gpJobStatus=function(jobInfo){
      console.log("Job Status");
      console.log(jobInfo);
  };

  var gpJobFailed=function(jobInfo){
      console.log("Job failed");
      console.log(jobInfo);
  };

In the above snippet the 'Layout_Template' is "A3_Landscape" and this works as i'm guessing the size of the map is predefined. I need to way to use the "Map_Only" template so that i can change the size of the map.


texts: gdal_translate error 1: output driver not found

texts:
I'm using gdalwarp from gisinternals msvc compiled binaries.

There should be a script which specifies the path variables. It was given inside the binaries. I modified in a small way such that it works from outside the directory and some other paths are also specified.

@echo off

set SDK_ROOT=%~dp0gdal\
set SDK_ROOT=%SDK_ROOT:\\=\%
set GNUPLOT_ROOT=%~dp0gnuplot\
set GSL_ROOT=%~dp0gsl2\

goto setenv
set PYTHON_ROOT=%~dp0python1\

@echo %SDK_ROOT%
@echo %GNUPLOT_ROOT%

if "%1" == "setenv" goto setenv

%comspec% /k "%SDK_ROOT%SDKShell.bat" setenv %1
goto exit

:setenv
@echo Setting environment for using the GDAL and MapServer tools.

if "%2"=="hideoci" goto hideoci

set ocipath=0
set _path="%PATH:;=" "%"
for %%p in (%_path%) do if not "%%~p"=="" if exist %%~p\oci.dll set ocipath=1

if "%ocipath%"=="0" goto hideoci
@echo WARNING: If you encounter problems with missing oci libraries then type:
@echo   SDKShell hideoci
goto setenv2

:hideoci
@echo Hiding the OCI plugin library.
if not exist "%SDK_ROOT%bin\gdal\plugins-optional" mkdir "%SDK_ROOT%bin\gdal\plugins-optional"
if exist "%SDK_ROOT%bin\gdal\plugins\ogr_OCI.dll" move "%SDK_ROOT%bin\gdal\plugins\ogr_OCI.dll" "%SDK_ROOT%bin\gdal\plugins-optional\ogr_OCI.dll"
if exist "%SDK_ROOT%bin\gdal\plugins\gdal_GEOR.dll" move "%SDK_ROOT%bin\gdal\plugins\gdal_GEOR.dll" "%SDK_ROOT%bin\gdal\plugins-optional\gdal_GEOR.dll"

:setenv2
@echo At set env 2
SET "PATH=%SDK_ROOT%bin;%SDK_ROOT%bin\gdal\python\osgeo;%SDK_ROOT%bin\gdal\python\;%SDK_ROOT%bin\proj6\apps;%SDK_ROOT%bin\gdal\apps;%SDK_ROOT%bin\ms\apps;%SDK_ROOT%bin\gdal\csharp;%SDK_ROOT%bin\ms\csharp;%SDK_ROOT%bin\curl;%GNUPLOT_ROOT%bin\;%GSL_ROOT%bin\;%SDK_ROOT%bin\gdal\plugins;%PATH%"
SET "GDAL_DATA=%SDK_ROOT%bin\gdal-data"
SET "GDAL_DRIVER_PATH=%SDK_ROOT%bin\gdal\plugins"
SET "PYTHONPATH=%SDK_ROOT%bin\gdal\python\osgeo;%SDK_ROOT%bin\gdal\python;%SDK_ROOT%bin\ms\python;%SDK_ROOT%bin\proj6\apps;%SDK_ROOT%bin\gdal\apps;%SDK_ROOT%bin\curl;%SDK_ROOT%bin;%SDK_ROOT%bin\ms\apps;%SDK_ROOT%bin\proj6\share"
SET "PROJ_LIB=%SDK_ROOT%bin\proj6\share;%SDK_ROOT%bin\proj6\apps"
SET "LD_LIB_PATH=%SDK_ROOT%gsl2\lib"

gdal_translate -gcp 0.0000 0.0000 80.207539 16.173769 -gcp 2683 0.0000 80.263502 16.173386 -gcp 0.0000 7906 80.206398 16.028650 -gcp 2683 7906 80.262312 16.028034 -gcp 1341 3953 80.234943 16.101114 HH_Amp.bin HH_Amp_geo.bin

pause

I specified the gdal_translate at the end of the script.

Now, I have specified the gdal_driver path in the script. To be sure, I also specified it in the path variable. But, where is this error coming from?


texts: Replace value in raster by another raster by ifelse (con) statement in QGIS

texts:
Previously I spent a lot of time to see answer on similar question  but neither work with my problem.
I have QGIS 3.4 and a DEM where I've to substitute more accurate height in cells that are identified by a "river" coming from rasterization of a poly line (... aligned and same cell width to DEM grid ..).
Outside the grid cell of river there are nodata cells, this is the problem I suppose, because ifelse function as 0/1 trick in QGIS doesn't work outside the river cells.
Also I'm not able to reclassify nodata in 0, neither of reclassification suggestion work.

When I was young ... I used the simple CON command in ArcView that worked in a very intuitive manner, to realiz this simple question "Where river raster exists pick height from this raster elsewhere maintain height of dem ..)
My administration believe to save cash to use open software but do not realize that I spent a lot of time to solve problem, as this, that as I know from previous experience should be intuitive.


texts: How do you prevent the data type in a CSV file from changing when you import the table from Excel into ArcGIS Pro?

texts:
I'm an ArcGIS newbie and I'm struggling to keep the data type from changing from numeric to text for my fields in a CSV table.

I've found that copying the table, pasting it into a new Excel document and then changing the fields manually there works, but I have to do that every single time I open the Excel table/import it to ArcGIS Pro. Also, I can't edit the field type in Pro and it's not letting me create a new field and set it equal to the original field I'm trying to edit.

I've included an example of the problem. What's also weird is the fact that only half of them get switched to text when I import it.


texts: Enable PostGIS extensions Windows local (non-admin) with Anaconda PostgreSQL

texts:
My organisation restricts admin privs on our machines, which means that installing OSGeo4Win would require a tech to come out (and as they're stretched pretty thin and I'm likely to need regular updates at least (but very likely several goes round as I mess this up) I'd love to find a solution that avoids them).

So I was interested to see that Anaconda bundles PostgreSQL (currently 11.2) as well as GDAL, PROJ4, GEOS etc. I can download the zipped PostGIS extensions, but am not sure where they ought to go.

Anaconda has put PostgreSQL executables (initdb, psql etc) in: C:\Users...\envs\geodjango\Library\bin\, along with a host of other binaries (many dlls.) that are not PostgreSQL.

The PostGIS zip says to drop the directory structure it contains into 'the postgres folder' & 'structure should mirror what you have.' So I copied and pasted the files and folders into "Library".  I am getting this error -

2019-08-16 14:14:39.495 AEST [10920] ERROR:  could not load library "C:/Users/...envs/geodjango/Library/lib/rtpostgis-2.5.dll": The specified module could not be found.

2019-08-16 14:14:39.495 AEST [10920] STATEMENT:  CREATE EXTENSION postgis;
ERROR:  could not load library "C:/Users/.../envs/geodjango/Library/lib/rtpostgis-2.5.dll": The specified module could not be found.

but the file is there:
c:...envs\geodjango\Library\lib>dir rt*

Directory of C:\Users...\envs\geodjango\Library\lib

15/08/2019  02:18 PM           740,864 rtpostgis-2.5.dll
               1 File(s)        740,864 bytes
               0 Dir(s)  164,782,927,872 bytes free

This is exactly the same file that was "missing" for this person:
Error in installing PostGIS extension through PostgreSQL 11


texts: Access all plugin from QGIS Modeler (Sextante)

texts:
I read that it is not possible to access to all of the installed QGIS plugins from the Graphical Modeler (QGIS 2.0.1).
It would be great, in defining the model in the Graphical Modeler, if the user could also add algorithms from all the QGIS plugins actually installed.

Is there a way to do this?

For example, I am interested in the use of the following plugins:

Add delimited text layer
Point Sampling Tool
Spatial Query


texts: How to get id of admin regions in Norway in OpenStreetMap?

texts:
I need to get the id of the relations that are administrative (boundary: administrative) that have admin_level equal to 7 in Norway.

Is there any way to get them?


texts: Are there torrents that contain large free spatial datasets?

texts: I've been having a lot of issues downloading entire datasets of SRTM or ASTER GDEM2 or similar free datasets, mainly because of poor services. Does anyone know of any torrents which contain such free datasets?


texts: Joining xls file with shapefile attribute through ModelBuilder?

texts:
I have folder contains xls files and I need to make join for each file with a shapefile attribute using ModelBuilder.

My trial is as follows:

and obtained "Precondition" option only.

How can I solve it?


texts: Openlayers create, edit and save/export to kml

texts:
I'm trying to create, edit and save a point vector as kml file in an Openlayers app. I have managed the create and edit by adding and empty vector layer and then adding the control

OpenLayers.Control.EditingToolbar("name of the layer")

But how could I save/export this afterwards in kml format? I just need this on the browser side a popup to let me save the file locally not on the server side.

Thanks
Aris


texts: what is the difference between the .mask and .updateMask

texts:
var l8 =ee.ImageCollection('LANDSAT/LC08/C01/T2_TOA')

var landsat2016 = l8.filterDate('2016-01-01', '2016-12-31');

var median = l8.filterDate('2016-01-01', '2016-12-31').median();

// Load or import the Hansen et al. forest change dataset.
var hansenImage = ee.Image('UMD/hansen/global_forest_change_2015');

// Select the land/water mask.
var datamask = hansenImage.select('datamask');

// Create a binary mask.
var mask = datamask.eq(1);

// Update the composite mask with the water mask.

var maskedComposite = median.mask(mask);

var maskedComposite2 = median.updateMask(mask);

So what the difference between maskedComposite and maskedComposite2?


texts: Combine named list of stars object into one new stars object with the names being a new dimension

texts:
I have a list of stars objects that I created in a loop. This little example  is more or less replicating my situation.
library(stars)
tif = system.file("tif/lc.tif", package = "stars")
x = read_stars(tif)
l = list(x,x)
names(l) = c("jan", "feb")

What I would like to have is one stars object with a temporal dimension that has two values: c("jan", "feb").
Now I know there is the generic c-function to combine stars objects and also the merge-function. But I did not manage to use the names of the list-elements as new dimension in the stars object. I don't event know if this is a good idea/approach. But in fact I quite often end up at this point (maybe I'm just doing it wrong all the time) and just wondered if it is easy to collapse a list of stars objects into one stars-object and use the names of the list as a new dimension.
EDIT
I did this:
s = c(do.call("c", l))
merge(s)

And this surprisingly creates the new dimension in the stars object, with as many values as elements in the list. But I'm not really sure if this is the best way...


texts:
For my use case, I'm specifically referring to HDF-EOS.

I'm trying to analyse ASTER data so am wondering which format I should utilise. I know GeoTIFF is easier and more compatible, yet HDF has everything needed and needed data can be efficiently sliced.

Although I've extensively read on the formats, I cannot seem to find a clear comparison for when to work with either. I also will be working with NetCDF in the future but I believe that's in the same bracket as HDF/EOS's strengths.


texts: Points and polylines change detection in QGIS (added, deleted, no change)

texts:
I'm relatively new to QGIS and am looking to find a simple plugin or method for doing the equivalent action as the "ChangeDetector" tool in ArcGIS. I have used Shapefiles previously with that tool. An example of how I used it:

I have two shapefiles
One is from 2016 the other from 2017
Comparing the 2016 to the 2017 shapefile, 50 new points have been added, 25 points removed and 500 with no change
The output from the ChangeDetector tool is a feature class called "New" with the 50 new points, another feature class called "Deleted" with the 25 deleted points and finally a 3rd feature class called "NoChange" with the 500 points

EDIT: when I think about it, I'm not really looking at spatial changes but rather how entries in the attribute table have changed, based on a unique ID number.


texts: Drawing Graphic at street layer/level only

texts:
I am new to ArcGIS JavaScript API and would like some help. I have a scenario where a user has to draw a graphic, mainly a polyline and some points on a map. To avoid confusion I want to limit the user to zoom in to street level prior to start drawing.

What I am trying to achieve here is when the user clicks on the locality he/she wants to draw on, on the click even the map zooms in to street layers.

I have been looking around but can't find a way to do it.


texts: Is there a way to clip a layer with a shape from another layer in `sf`?

texts:
I would like to get the greens polygons that intersects (or overlapsed) the red outlined polygon (see image below).

I tried sf::st_intersection() but the polygon in yellow (yes that just one polygon) was not taken in account.

In QGIS, the clip (cut) tool did exactly what I wanted. Is there a tool similar in R to cut a layer using another layer ?

Tried to use sf::st_overlaps() but the returned matrix was biggest than my feature table. sf::st_crops() does not fit either as I want to cut in a specific shape.


texts: Python gdal.translate doesn't write the output file

texts:
I am very new with gdal and python in general, therefore pardon my simple question.
I have a geotiff (single band, black and white) that I want to change from grayscale to RGBA (becasue I want to use gdal.DEMProcessingOptions to add color-relief).
I am using the following script but I don't get any output file in the folder. Do I have to specify somewhere that I want to write the output?
gray = gdal.Open("input.tif")
rgba= "output.tif"
colour=gdal.TranslateOptions(rgbExpand="rgba")
gdal.Translate(rgba, gray, options=colour)