texts: ArcGIS10 VBA requires library references

texts:
I just upgraded from ArcGIS 9.2 to 10.0, and I'm trying to create some tools using VBA. It doesn't recognize any of the ArcObjects Classes, Interfaces, or Constants unless the reference includes its library name. For example:

Dim pMxDoc as IMxDocument
Set pMxDoc = ThisDocument

doesn't work, but

Dim pMxDoc As esriArcMapUI.IMxDocument
Set pMxDoc = Project.ThisDocument

does.

Does anyone know how to fix this? (FYI I have no admin capabilities on this computer.)


texts:
I'm trying to follow a series of tutorials:

visualizing-terrains-with-cesium-ii

visualizing-terrains-with-cesium

Some other useful posts here:

Unable to serve Cesium Terrain files

error-in-loading-terrain-files-terrain

Basically, I want to use my own terrain files within Cesium.

I'm stack in configuring the terrain server aspect,

I have a Node.js server running Cesium Apps in port 8080:
http://localhost:8080/Apps/own_terrain.html

Path for the html file is:

C:/Users/username/Downloads/software/Cesium-1.53/Apps/own_terrain.html

The terrain files have been generated by using: cesium-terrain-builder-docker

Then, I want to run my terrain tiles using cesium-terrain-server in docker.

This is how I built it:

docker run -p 9000:8080 -v C:/Users/username/Downloads/software/Cesium-1.53/Apps/terrain:/data/tilesets/terrain \
geodata/cesium-terrain-server

So, Cesium running in 8080 and terrain server in 9000,

However when I try to access the terrain tiles, i.e:

http://localhost:9000/tilesets/tiles/0/0/0.terrain

Nothing happens, Chrome says:

localhost didn’t send any data. ERR_EMPTY_RESPONSE

And the same when integrating that into the own_terrain.html:

Cesium.js:445 GET http://localhost:9000/tilesets/tiles/layer.json net::ERR_EMPTY_RESPONSE

Definitely the terrain server is not running. What am I missing in these tutorials?


texts:
I made an exe file from standalone application. However I want to enter the path variable on cmd when running the exe.

My standalone application running *.shp files. But I want to enter path of the *.shp files at the beginning, when I run the exe on cmd.

For instance my exe file name is myMap.exe and I want to do that:

C:\Users\Mustafa Uçar>cd C:\Users\Mustafa Uçar\Desktop\MyMap\build\exe.win-amd64-3.6

C:\Users\Mustafa Uçar\Desktop\MyMap\build\exe.win-amd64-3.6>myMap.exe C:\Users\Mustafa Uçar\Desktop\Tutorial\Mississippi\mississippi.shp

PyQt5 - Python 3.6 - Qgis 2.99

Here is my code (I give path of shapefile but I do not write it in script. I want to write it on cmd. Path variable must be changeable):

import sys
from qgis.core import *
from qgis.gui import *
from PyQt5.QtCore import *
from PyQt5 import QtGui

app = QgsApplication([], True)
path = "C:/OSGeo4W64/apps/qgis-dev"
app.setPrefixPath(path, True)
app.initQgis()

canvas = QgsMapCanvas()
title = "PyQGIS Standalone Application Example"
canvas.setWindowTitle(title)
canvas.setCanvasColor(Qt.white)

filepath = "C:/Users/Mustafa Uçar/Desktop/Tutorial/Mississippi/mississippi.shp"

layer = QgsVectorLayer(filepath, "testlayer", "ogr")
QgsProject.instance().addMapLayer(layer)
canvas.setExtent(layer.extent())
canvas.setLayers([layer])
canvas.zoomToFullExtent()

canvas.show()
sys.exit(app.exec_())


texts:
I digitized some maps with isolines every 4m. I created a DEM and the Hillshade. But the hillshade has big pixels. I use ArcGIS 10.3.

Is there something to be look out for that I might have done wrong?


texts: Save lappy results without overwritting - R

texts:
Here is my situation,

I have 10 multi-spectral images. For each one I am calculating three different vegetation indices. Those 3 raster objects share resolution, projection, extent and dimension (they are produced with raster calculations).

I am using this code to list them (not stack since I will use lapply later)

VI<-as.list(Index1,Index2,Index3)

The next step in my methodology is to upload a shp representing the field boundaries to derive the mean, min and max value inside each polygon FOR EACH raster of the list.

The code I am using for this is:

shp<-readOGR("myshapefile.shp")

Statistics<-function (VI) {

    meancalc<-extract(VI, shp, fun=mean)
    mincalc<-extract(VI, shp, fun=min)
    maxcalc<-extract(VI, shp, fun=max)

    shp@data<-data.frame(shp@data, mean=meancalc)
    shp@data<-data.frame(shp@data, min=meancalc)
    shp@data<-data.frame(shp@data, max=meancalc)

   writeOGR(shp, getwd(), layer="mylayer", driver="ESRI Shapefile")
}

I run the function using
    output<-lapply(VI, FUN=Statistics)

The output is the shp with three new attributes (mean, min and max). However, it is obviously overwriting the calculations each time lapply takes a new raster as input for the function.

What I want is the shapefile with individual attributes for each raster input. Something like: Index1mean, Index1min, Index1max + Index2mean, Index2min, Index2max + Index3mean, Index3min, Index3max

I know a solution would be to link the output name of the shp to the input raster name so that three different shp are created for each iteration. However, I would prefer to keep everything in the same shp

Any idea how to implement this on the code?


texts: Automated merge of polygons up to a limit

texts:
I have a polygon feature class with a field, let's say it's population. I need to find an automated way to merge these polygons up to a given population number.

E.g. I have polygons with population from 0 to 211. And I need to create districts with about 400 habitants. (Between 380 and 420.)

I work with ArcGIS 10.


texts: Using ArcPy UpdateCursor to Update a field

texts:
In my data I have one field (L1_SOILTYP) and a new field I created (SoilClass). I am trying to use update cursor to assign classification values (1-17) based on the values found in the L1_SOILTYP field. For example... 'Clay' and 'CLAY' should be assigned '1'. I am using ArcGIS Pro for this.
Following an example on this page... https://pro.arcgis.com/en/pro-app/latest/arcpy/data-access/updatecursor-class.htm ... this is the code I have come up with...
import arcpy
fc = 'D:\GIS\PollinatorProject\Pollinator\Pollinator.gdb\soil_samples'
fields = ['L1_SOILTYP', 'SoilClass2']

with arcpy.da.UpdateCursor(fc,fields) as cursor:
    for row in cursor:
        if(row[0] == 'CLAY,DUMP'or row[0] == 'CLAY(WEATHERED SHALE)'or row[0] == 'FILL,CLAY'):
            row[1]= '1'
        elif(row[0] == 'SAND' or 'TOO SANDY' or row[0] == 'Sand, coarse' or row[0] == 'COARSE SAND'):
             row[1]= '2'
        elif(row[0] =='Sandy Clay' or row[0] == 'SAND & CLAY' or row[0] == 'SAND SOME CLAY'):
            row[1] = '3'

        cursor.updateRow(row)

Using this code the clay values are all correctly assigned the 1 value, but all the other soil types are assigned to 2. This is a shortened version of what I have, but each time I run the code, all values are assigned '1'


texts: Raster layers do not appear in qgis2web layer list

texts: Using qgis2web I created a Leaflet map that has a vector layer, a raster layer and OSM background. I included in the map the layer list ("Add layer list" option in qgis2web). The vector layer appears in the layer list and I can turn it on/off. However the raster layer does not appear in the list. Is that the intended behavior? Using OL3 I can turn rasters on/off.


texts: Labeling a point layer in QGIS: What is the "around point" option doing

texts:
When you want to label a point layer in QGIS and go to "Position" there are the three options shown in the image (sry its not in english):

Hovering on the "around point" option gives me the information

  label will be positioned on a circle around the point object

I don't understand what that means. Anyway we just have one label right? Thus it will be positioned on a circle around the object but I won't get many labels around the object. What could be usecases for this option?

I'm using QGIS 3.4.13


texts:
I am trying to create a base map layer of a district in Leipzig, Germany. It should consist of building shapes and addresses and house numbers so my team can add the data we are currently collecting. They have certain values in a excel sheet with referenced addresses.

I downloaded OSM data from Geofabrik. However, this data also doesn't have all addresses.

This is the data I have in QGIS. Houses with building numbers are colored.
I just don't understand how the online OSM map has all the adresses and the one I downloaded from Geofabrik don't.

The same view in online OSM

Could you point me to a direction please? I have tried several approaches such as downloading the whole GIS data for Saxony from Geofabrik, using overpass-turbo to query addresses but it seems that the map sometimes stores the address in nodes and sometimes in ways. How could I combine them?


texts: Looking for advice on getting GIS data into existing AutoCAD site survey

texts:
I've got a variety of GIS data from Albemarle County in ShapeFile and DWG format. I'm looking to get the topographic layers and property boundaries together with my existing site survey created by a local surveyor.

What's the best practice for doing this? Is it "just XREF it all into one master file" and scale/move around or is there a better way to make this data normalized? I'm in one of the corners of the grid squares that the county has created, so I'd probably need to have 4 of them combined.

Edit: I have AutoCAD 2013. I have very little true GIS experience, but I have friends who have offered to help with QGIS.


texts: Setting up account to start using Collector for ArcGIS?

texts:
I'm working on a presentation topic around mobile-based data collection solutions. I've used ODK and Fulcrum in production and absolutely love them. However, I'm going to be presenting to an Esri crowd so I want to at least add Collector for ArcGIS to the pool of candidates for comparison.

How do you get this working?

It's difficult to "test"  this product out. I've had an Esri global account now for about 10 years.  I can't use that ArcGIS.com account in collector because it is not an "organizational account" (which I don't understand because it is the login tied to our corporate licenses). I created a new account using the "30-day Free Trial" found at http://www.esri.com/software/arcgis/arcgisonline.

I created a new map, added an editable layer to the map, saved the map, shared the map publicly, and to a fake group I made, I even added a dummy email (user account) as a member of the group/organization, and yet when I go into Collector and login with either my admin email, or the dummy email, no maps show to be available.

I've been poking around online and within the Esri forums this morning for about 3 hours trying to troubleshoot the issue. My first thought was that I couldn't use Collector with a "trial" ArcGIS Online account, but most everything I read from Esri reps stated otherwise. Perhaps that has changed.

Does anyone have any insight to provide on this issue?


texts: Create a marker in the bottom right corner of the map

texts:
I try to set a marker when on the bottom right of my map no matter which zoom level i have, similar to a control where i can set the position to 'bottomright'. What i want to achieve is to show a circle everytime a marker is dragged so i am able to drop the marker in that circle. What i tried so far is this:
private createDropZone() {
    const mapBounds = this.map.getBounds();

    const x = mapBounds.getSouthEast();

    const lat = x.lat - 0.02;
    const long = x.lng - 0.02;

    this.dropZone = new L.CircleMarker(new L.LatLng(lat, long), {
      color: '#f27e00',
      fillColor: '#f27e00',
      fillOpacity: 0.3,
      radius: 50,
    });

    this.dropZone.addTo(this.map);
  }

So as i you can i am subtract 0.02 from the south east bound of the map but of course this will only work on one zoom level. If its not possible to draw it on the same position in the right bottom corner it would help if someone give me hint how i can draw a marker relativ to the given marker.


texts:
I have a geological style file Biblioteca_CPRM_v2.

But when I open it it appears like this

Should look like with this

I guess it may be something in style fonts. Any idea?


texts: Arcpy InsertCursor inserting geometry into Shape field

texts:
I'm trying to write a script that reads geometry for each row in feature class, then create a new feature class and insert into it that geometry. It has to work with ArcGIS 10.0 so I have to use arcpy.InsertCursor. But I have a problem, because in the new feature class there is the same number of rows as in input but every row has the same geometry, from last row of previous feature. In 10.2, even using arcpy.InsertCursor, everything is working fine.

import arcpy
arcpy.env.workspace
infile=arcpy.GetParameterAsText(0)
outfc="przekroje_oblicz_hydroo.shp"
outpath="C:\\Users\\Bartek\\Desktop\\temp"
shapeName = arcpy.Describe(infile).shapeFieldName

linie=[]
km=[]
ID_HYDRO_R="12334"
#ver=arcpy.GetInstallInfo()['Version']
ver="10.0"

if ver=="10.0":
    cursor = arcpy.SearchCursor(infile)
    for row in cursor:
        linie.append(row.getValue(shapeName))
        km.append(row.getValue("Km"))

else:
    for row in arcpy.da.SearchCursor(infile, ["SHAPE@","Km"]):
        linie.append(row[0])
        km.append(row[1])

PUWG92 = "PROJCS['PUWG_1992',GEOGCS['GCS_ETRS_1989',DATUM['D_ETRS_1989',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['false_easting',500000.0],PARAMETER['false_northing',-5300000.0],PARAMETER['central_meridian',19.0],PARAMETER['scale_factor',0.9993],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]]"
fc=arcpy.CreateFeatureclass_management(outpath, outfc, "POLYLINE", "", "", "ENABLED", PUWG92)
#fc=arcpy.CreateFeatureclass_management(outpath, outfc, "POINT", "", "", "ENABLED", PUWG92)
arcpy.AddField_management(fc, "ID_HYDRO_R", "STRING", "", "", "22", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(fc, "KM_PKT", "DOUBLE", "6", "3", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(fc, "ZW_10", "DOUBLE", "7", "2", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(fc, "ZW_1", "DOUBLE", "7", "2", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(fc, "ZW_02", "DOUBLE", "7", "2", "", "", "NULLABLE", "NON_REQUIRED", "")
arcpy.AddField_management(fc, "ZW_WZ", "DOUBLE", "7", "2", "", "", "NULLABLE", "NON_REQUIRED", "")

if ver=="10.0":
    k=0
    cur = arcpy.InsertCursor(fc)
    for i in linie:
        feat = cur.newRow()
        feat.Shape = i
        feat.ID_HYDRO_R=ID_HYDRO_R
        feat.KM_PKT=((float(km[k])/1000))
        feat.ZW_10=11
        feat.ZW_1=12
        feat.ZW_02=12
        feat.ZW_WZ=-7777
        cur.insertRow(feat)
        k+=1

else:
    cursor = arcpy.da.InsertCursor(fc, ["SHAPE@","KM_PKT"])
    k=0
    for i in linie:
        kma=(float(km[k])/1000)
        #arcpy.AddMessage
        cursor.insertRow([i,kma])
        k+=1

del cursor, cur

EDIT:

I managed to rewrtie geometry by modyfying searchcursor to read vertex into 'arcpy.Array', then create 'arcpy.Polyline' and add it to list linie. But it is much slower solution.

cursor = arcpy.SearchCursor(infile)
for row in cursor:
    feat=row.getValue(shapeName)
    k=0
    punkty_lista=arcpy.Array()
    for part in feat:
        for pnt in feat.getPart(k):
            punkty_lista.add(pnt)
    linia=arcpy.Polyline(punkty_lista)
    linie.append(linia)
    punkty_lista.removeAll()
    k=+1


texts: How to modify or create a wider UTM zone?

texts:
am working on a project to do country-wide mapping using projected coordinate system such UTM. However, the country (Zambia) is covered by three UTM zones, namely zone 34,35 and 36.I have read literature on possibility of widening UTM zones and the resulting worsening of errors.
My concern is on how to modify UTM zones to cover as much as 12 degrees of longitude in one zone. The central meridian will of course be moved to the center.


texts: How to Drop Indexes on a Table in ArcCatalog

texts:
I am using an Oracle 11g database . The code below doesnt allow me to drop indexes on a table.

import arcpy

# Remove indexes and statistics for gis feature class delta tables
COND=r'E:/GISScheduledTasks/DMZGeodatabaseMaintenance/sde/arcfm_synctest.sde/ARCFM.CONDUCTORINFO'
arcpy.RemoveIndex_management(COND,"G190PRIOHCONDUCT;G190PRIOHELECTRI;G190PRIUGCONDUCT;G190PRIUGELECTRI")
print('complete')

It gives the attached error.Could someone advise on python coding for dropping indexes on a table within an enterprise database?


texts:
On https://gis.stackexchange.com/a/19979/19654 i found a solution to add hover / mouseovers to my map:

       // add code to create tooltip/popup
       popup = new OpenLayers.Popup.FramedCloud(
          "",
          feature.geometry.getBounds().getCenterLonLat(),
          null, //new OpenLayers.Size(320,20),
          "<div>" + feature.attributes['name'] + "</div>",
          null,
          true,
          null);

       feature.popup = popup;

Doing this i found the popup don't have a nice / clean layout, see:

I tried with a fixed size too (new OpenLayers.Size(320,20)) which seems to have the same effect.

So, my question would be if someone has a cleaner/easier/nicer way to improve this? Thanks in advance for any hint on this.


texts: PyQGIS - Adding Fields to a Feature Layer

texts:
I'm new to PyQGIS and am having difficulty adding fields to a feature.  I want to plot footprints and assign attributes to each of the footprints.  Here is my script.

import csv, ast
from PyQt4.QtCore import QFileInfo, QVariant

fileName = 'C:\\directoryname\\temp.csv'
layer = QgsVectorLayer('Polygon?crs=EPSG:4326', 'poly', 'memory')
pr = layer.dataProvider()
attr = pr.addAttributes([QgsField('Object_ID', QVariant.String)])
poly = QgsFeature()
with open(fileName) as f:
    rows = csv.DictReader(f)
    for row in rows:
        shape = ast.literal_eval(row['footprint'])
        points = []
        for i in range(0,len(shape['coordinates'][0])):
            points.append(QgsPoint(shape['coordinates'][0][i][0],(shape['coordinates'][0][i][1])))
        poly.setGeometry(QgsGeometry.fromPolygon([points]))
        poly.setAttribute('Object_ID', row['object_id'])
        pr.addFeatures([poly])
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([layer])
print "Completed!!!"

I receive the error

  File "c:/directoryname/tmp9urit7.py", line 19, in

      poly.setAttribute('Object_ID', row['object_id']) KeyError: 'Object_ID'

When I add the line caps_string = layer.dataProvider().capabilitiesString(), then printing caps_string yields that I have full privileges to add, modify, and delete features and attributes.

Also, when I print attr, the result is True.

If I comment out the line poly.setAttribute('Object_ID', row['object_id']), the script plots the footprints and completes without error, but I have no attributes for any of the footprints.

What am I doing wrong?


texts: finding the definition of missing srid in postgis

texts:
I am converting a spatial Oracle database to POSTGIS, and one of the problems I encounter is that we use SRID in Oracle that do not exist in POSTGIS.

I am talking about the following SRID:

90112 : dutch RD
327680: Belgian Lambert 72
8307: WGS84

The problem is that I do not know how to translate these SRID to POSTGIS, since they do not exist in spatial_ref_sys.

If I visit the spatial reference I do not find any reference to the SRID at all. I was able to find

31300: Belgium Lambert 1972
4326: WGS84

Are these SRID comparable? Where can I find definitions to add SRID if they are not defined on spatialreference.org? Are the Oracle SRID's standard or not?

A lot of questions, I hope you can shed some light here. Thanks :)


texts:
I have polygons from coordinates in (python shapely) that looks like this

POLYGON ((24.8085317 46.8512821, 24.7986952 46.8574619, 24.8088238 46.8664741, 24.8155239 46.8576335, 24.8085317 46.8512821))

I would like to calculate the area of this polygon in km^2. What would be the best way to do this in Python?


texts:
I am trying to save my openlayers map on the server, following these instructions, the sections Client-side: Tile collection and Server side: PHP:

http://trac.osgeo.org/openlayers/wiki/TileStitchingPrinting

On the server side, the script creates lots of tiles-files, with names like 6910a8c32e7e8e4e81270bba8a106231.img, and, finally, an jpg file , 3696d338f25f6a64ea6f90a78b3276b9.jpg, that is supposed to be the image of the maps. However, the content of this file is only black, whatever the content of the map might have been.

Any ideas?


texts:
I am trying to find a good way to serve ortophoto with geoserver with fast rendering.
I have 20GB of data, 4000 jpeg tiles at 6-15mb each.
I have published them as a mosaic but rendering them is very slow as a wms and I have not yet been able to preview tile layer (server never respond - Handler dispatch failed; nested exception is java.lang.OutOfMemoryError: GC overhead limit exceeded). Are the files to many/too large? Are there any recommendations how to serve a image mosaic?
In the past I have been creating an imagePyramid with 4 levels first, which I then have cached, but that feel a bit redundant since I am mostly interested in the cache. Creating the imagepyramid also reduced the image quality too much at high zoom level. I was hoping that publishing the image mosaic would keep the image quality better for the cache at higher zoom levels.
What is the best approach to take to retain the image quality on higher zoom level and at the same time be able to serve the tiles quickly? I was hoping I could use a image mosaic so i can easily replace individual tiles easily when i have updated data.


texts: Arcpy Mosaicing Rasters - Output seems to be wrong

texts:
I am using the Arcpy-package for Python 2.7 (with Arcgis 10.1). I tried to write a script that is mosaicing all rasters in my directory, where they overlap it should calculate the mean for the rasters.

This is my script:

import arcpy

arcpy.env.workspace = "D:/RSIV/Masterarbeit/GRD/Norway/Neu/SentinelSat"

# Import all raster with .TIF in the environment.
rasters = arcpy.ListRasters("*", "TIF")
arcpy.Mosaic_management (rasters[1:],rasters[0] , "MEAN","FIRST","0", "0", "", "0", "")

It works fine, at least its producing the result and overwriting the first raster, but the values are simply wrong for some parts of the image.
The image below shows the identify-tool from ArcMap 10.1, as you can see there is no data in 2 of 3 bands, but the Mean is a little bit of the expected value.

If someone can tell me where I have to update my script/what is going on with the data I would be super thankful.
Note: I do not have to use Arcpy for this, if there are other packages I am willing to use them, even for python3.X.
(Right now I am thinking about Rasterio but I did not work with it before, will give it a read tomorrow)


texts: QGIS polygon border width option has disappeared

texts:
After upgrading to QGIS 2.16, I can no longer see any option to change the border width of polygons in the 'style' option in properties. I want a transparent polygon with a thick black border.
I choose 'single symbol' then 'outline: simple line'.
But I can't see any options to change the colour or width of the border, which was always there in previous QGIS versions.
Instead, it says 'this layer doesn't have any editable properties'.
Am I missing something obvious?


texts: Retrieving geoprocessing results in ArcGIS JavaScript API

texts:
I'm executing a geoprocessing script in an ArcGIS JavaScript API map viewer, and it works until I attempt to read the job's output parameter.

My script:

function extractModel(results) {
    var params = { "ModelName": ModelName }; // variable defined earlier

    // run geoprocessing service to export GDB, download results
    var gp = new Geoprocessor("http://myserver/arcgis/rest/services/Test4Custom/GPServer/Custom");
    gp.setUpdateDelay(5000); // check status every 5 seconds
    gp.submitJob(params, statusDone, statusCallback, errorBack);
}

function errorBack(jobInfo) {
    alert.window("Error encountered in geoprocessing script.")
    console.log("Status: " + gpStatus);
}
function statusCallback(jobInfo) {
    console.log("Status: " + gpStatus + " -- Continuing...");
}
function statusDone(jobInfo) {
    console.log("geoprocessing completed");
    console.log("Status: " + gpStatus);
    console.log("job id " + jobInfo.jobId);
    gp.getResultData(jobInfo.jobId, "Output_File", downloadResult, errorBack);
}

function downloadResult(result) {
    console.log("displaying result");
    console.log(result.value);
    console.log(result.dataType);
}

The geoprocessing task executes correctly, and the viewer script is running up until the gp.getResultData line. It's not getting into the downloadResult() function, and I'm getting this error in the console:

  TypeError: c is not a function(…) "TypeError: c is not a function

Does this indicate a problem with my syntax in calling downloadResult(), or something else?

Note: The output parameter is indeed Output_File, and has a string in it when executed. Example:

{
 "paramName": "Output_File",
 "dataType": "GPString",
 "value": "c:\\arcgisserver\\directories\\arcgisjobs\\test4custom_gpserver\\je47fc8d361064c9e86237f397a873118\\scratch\\Data_SALLEY_20151215.gdb"
}


texts: In eCognition developer is it possible to save samples that I have selected for nearest neighbor classification with the workspace or project? Every time I close and reopen my project I have to re select NN samples.


texts: R process: Zonal Stat Raster has No Data Values vs Arc output

texts:
I am trying to re-create my ArcGIS Zonal Stat model with R
All my files will be 30m pixels; 18000 x 18000.
I have used R two different ways to create the output and both create NoData values in the bottom portion of the output raster.
My Arc GIS Zonal Stat Max output has values of zero & one.
My R zonal Stat Max is mostly zero & one except for a portion of the image; ~~ 86 pixels from the bottom up, across the entire bottom, primarily come out as zero & No Data.
I have included a picture of my output rasters showing the one section within those 86 pixels that actually write values of zero. I cannot just flip my NoData values to zero as there are other sections that need a value of one that are also getting set to NoData.

Python 2.7 with rpy2 created the clump file.
To create the Zonal File, I needed to write an R script file then CALL R within python.
I used the link below to write my R script.

http://www.guru-gis.net/efficient-zonal-statistics-using-r-and-gdal/

system.time(q <- reclassify(z, as.matrix(Zstat)[,1:2]))

Where z is my zonal file & Zstat is my zonal data frame.

After running line by line in R, I figured out the line of code above writing of the output file was taking the longest. 12 hours to complete or create the (q variable) on a 18000 x 18000 30m file PLUS NoData values in the bottom ~ 86 pixels of the output.
I then changed the code to create a matrix from the data frame with code

I found here:
    https://stackoverflow.com/questions/52073767/joining-tables-to-a-raster-attribute-table-to-create-new-raster-with-the-reclass

system.time(reclass <- matrix(Zstat[,"vals"][as.vector(z)+1]))

This now takes 20 min & another 20 min to create (q - below)

system.time(q <- raster(reclass, template = z))

When calling both of these in the python script it is about an hour vs 40 min within Python, line by line in R is less than a minute a piece.

The main problem other than the time it takes is that my output raster,  as mentioned, for the last 86 pixels contains NoData values.

which(is.na(Zstat))
results: integer(0)

HOWEVER, after the matrix as vector.

((reclass <- matrix(Zstat[,"vals"][as.vector(z)+1])))

is run.
I now have NoData values where I should NOT……

head(which(is.na(reclass_2)))

>[1] 322595978 322596008 322596041 322596042 322596118 322596237

> length(reclass_2)
[1] 324036001
> na <- which(is.na(reclass_2))
> length(na)
[1] 32901

So what is causing 32901 pixels to be NA..?

Images:

First: right most green dot, in grey NoData area.

Second: left most green dot, in brown value zero area

third: dark dot in green area, green is NoData and should be value 1

Files mentioned:

Final_Zonal_Stat_file: R zonal output raster

ARCgis_Zonal output: created with ArcGIS model

Clump/Zones File: created in python with rpy2

water information: python gdal recoded raster file

UPDATE:

system.time(reclass <- matrix(Zstat[,"vals"][as.vector(z)+1]))

The clump file is missing sequential values for some reason, even though R was set to gaps=FALSE..... so my output for Zstat looks like this in the problem area:

         z    vals
1872968 1872967 0
1872969 1872968 0
1872970 1872969 0
1872971 1884697 0
1872972 1884698 0
1872973 1884699 0

The jump in numbers in row 1872971 is the problem. All the second column from this point on are not written to the final output raster. or reclass variable. and result in NoData values in my final matrix/raster.


texts:
I am trying to add an HTML table to the bottom left corner of my print composer. The table is populating with data based on the current atlas feature. The problem is that the table has a different number of rows based on the atlas feature. The way the html table renders in composer currently it draws from the top left corner down. This means that if the box is centered for two rows and an atlas feature has a different number it is no longer lined up properly. Is there a way to have the HTML table draw from the bottom left corner up? This way the table would still be placed properly on the page no matter the number of rows.

The only solution I have found so far is to have the whole composer draw upside down and then rotate the map back to the correct orientation after the export.


texts: How do I increase the font size of contents in Map in Leaflet.js?

texts:
I am using Leaflet JS for displaying my map.

Following is my code. The map gets displayed. But the size of the city names are not large. Each time I zoom in, the font size is very less. It is not readable.

var mymap = L.map('maps1').setView([19.8762, 75.3433], 7);
    L.tileLayer('https://api.mapbox.com/styles/v1/{id}/tiles/{z}/{x}/{y}?access_token={access_token}', {
    attribution: '&copy; OpenStreetMap contributors',
        id: 'mapbox/streets-v11',
        accessToken: 'your.mapbox.access.token'
}).addTo(mymap);
L.marker([19.8762, 75.3433]).addTo(mymap).bindPopup('Aurangabad');


texts: Download OSM data with QuickOSM and QGIS Modeler

texts:
I've been using QuickOSM plugin for some time and I would like to create some model that allows me to automate OSM data dowload in order to make sure that I am always working with the most updated information.

My intention was to perform an overpass query on bbox defined by a custom name to look at nominatim and then save the output into a vector file and load it into QGIS in order to work with it. So, with that in mind I created the following very basic model (it's the first one I ever create -a part from a basic tutorial example):

(check code here)

Unfortunately after creating this basic model not only I don't get any results but I get an error "Query overpass API with a string Bad request OverpassAPI". Additionally, I cannot document my model by adding a title, description... so I am asuming I did something wrong. Additionally I guess that my model will only save a file containing points but will not load it to QGIS as I would like to.

Can anyone check my model and/or guide me to achieve such a basic behaviour?


texts: Creating Atlas dynamic text from multiple layers in QGIS

texts:
I am creating an atlas using QGIS v3.12 but would ideally like to include a dynamic text box that takes values from multiple layers contained on each atlas feature page.
The map is made up of two layers ('zones' & 'networks')
'zones' used to generate the atlas, but I'd like to include some of the attributes from 'networks' including 'total_properties'.
I've tried using:
aggregate(
layer:='networks',
aggregate:='concatenate',
expression:="total_properties",
filter:=intersects($geometry,@atlas_geometry),
concatenator:=','
)

But it returns the error "Eval Error: Could not calculate aggregate for: total_properties".
It only seems to work if the attribute name is common with the 'zone' layer (they are different geometry types from different sources but have some common attributes). I even created an attribute table for 'networks' and copied the attribute names from the expression builder to ensure they're spelt correctly.


texts: Exporting gdb tables to excel in folder containing gdb?

texts:
So i have several folders which each contain a gdb.

Now, I export the attribute tables of polygon feature classes using the following code (probably could use some optimization):

folder = arcpy.GetParameterAsText(0)

for dirpath, dirnames, filenames in arcpy.da.Walk(folder, datatype="FeatureClass",type="Polygon"):
    for filename in filenames:
        fc = os.path.join(dirpath, filename)
        table = dirpath + os.sep + filename + "_Table"
        result = arcpy.GetCount_management(fc)
        if int(result.getOutput(0)) > 0:
            print "Creating " + table
            arcpy.CopyRows_management(fc, table)

Output looks like:

-folder
--subfolder1
---*.gdb <--contains multiple tables
--subfolder2
---*.gdb <--contains multiple tables

I would like to simply output an excel version of the table at the same time, but i am unsure how to parse out the folder containing the .gdb to use as the path for the excel file.

Outout should then look like:

-folder
--subfolder1
---*.gdb <--contains multiple tables
---Excel1
---Excel2
---etc.
--subfolder2
---*.gdb <--contains multiple tables
---Excel1
---Excel2
---etc.

How would I go about getting that parameter to output the excel file in the correct location?


texts: Why does Clip function (maintain clipping extent) not clip according to polygon feature?

texts:
In arcGIS 10, I have several raster maps, that I would like to clip according to a polygon feature (top map). I have tried the Clip function under Raster Processing<-Raster<-Data Management Tools, and ticked off the Maintain Clipping Extent box. I am not allowed to tick off the Use Input Features for Clipping Geometry box.
But I end up with a map of the same extent but different shape than my polygon (bottom map). I really need the map to have the same extent and size of the polygon. I have tried un-ticking the Maintain clipping Extent box, but it makes no difference with regards to the result.
I have added a picture of what I entered in the Clip box, and what the process looked like after the maps.


texts: How to use ArcGIS WPF Runtime Deployment Tool via Command Line?

texts:
I'd like to be able to use the WPF Runtime (10.1.1) Deployment tool via the command line. If I run DeploymentTool.exe /? on the command line, I get the dialog in the attached image. The /c option indicates support for using a config file. I cannot find documentation on the config file format. Can someone point me in the right direction?

Apparently when you create a configuration in the tool, it is stored in the registry at HKEY_CURRENT_USER\Software\ESRI\ArcGIS Runtime SDK\Settings. If any additional information is stored elsewhere, I can't find it.


texts:
How could i create e mapping applications that displays a custom map from this URL:
(Capabilties with a mapview)
http://openwms.statkart.no/skwms1/wms.sjo_hovedkart2?request=GetCapabilities&Service=WMS

inside c# with WPF controls.

It would be nice if there was some full feature map libraries so i could add pinpoints to certain Long,lats e.e (Like Openlayers for Javascript)

Google / bing maps API cannot be used as you are not allowed to replace their maps.


texts: Combining python and JS (mapbox) for interactive rendering

texts: I plan on visualising a huge dataset (probably around 20-30GB of data) using Mapbox in order to have a very nice map visualization. However I would like this visualization to be interactive and I am pretty sure that both JS and my browser will not really like such an amount of data, so I wanted to have my data filtering/processing operations conducted with python and the rendering done dynamically with JS and Mapbox. I have seen that there is Python API for Mapbox, but all of the examples I have seen so far of it simply create a geojson dataset once and a html page, but it seems that the rendering code in JS is not dynamically generated. Would you guys know if there is any way to achieve that?


texts: QGIS 3x label appearance at different zoom level

texts:
I have a problem.
In my QGIS project, all the labels appear the same despite the zoom level. What looks good when the zoom is big, it's not as great when zooming out.
I found a way to get it disappeared,
Label display at different zoom level of single shapefile?
but I need it smaller only. Is there any way to set the different font sizes for different zoom levels?
I've already the expression based on one of the table value:
  CASE
  WHEN "Visible" = 1 then color_rgb (0, 0, 0)
  End

How can I integrate this expression with the answer below?


texts: List of elements in join layer without a match

texts:
When I set up a table join in one layer A, how do I find all elements from layer B that didn't find their way into the joined table?
This is much like doing a left outer join and trying to find all elements from the right table (B above) that didn't go into the join.
I suppose I should be able to set up a join in layer B and use layer A as the lookup table. Then, I could just mark all elements with NULL values in the joined fields, but somehow, QGIS won't let me set up a join in layer B for layer A. Which is really strange.

Layer A has 5000+ features.
Layer B has 280 features.
When I inspect the attribute table for layer A and sort layer A by an identifier of layer B, I get 268 joined elements. OK, I have 12 elements in layer B without a match. Of course, since they don't appear in the joined attribute table, I would need to look for them in layer B. But in layer B, of course, they are unmarked.
Setting up a join in layer B properties for layer A doesn't work (the join settings window remains blank after setting up the join). Well, then I tried "join by attribute" and did not skip the faulty output, which I saved as separate layer. This layer should have 12 elements, but only has 10. I must be doing something wrong. How would you find the unjoinable features in layer B?


texts: VB.net Converting UTM to Lat/Lon (State Plane with WKID) - ArcGis

texts:
I have looked at many libraries. Many simply don't work, others are very poorly documented.

I have a very long list of  X/Y coordinates (US FT) and an associated WKID.  I am trying to convert that to Lat/Lon format for mapping.

Can anyone suggest routes for completing this process in vb.net?

(the closest I've come to success is finding the esri.arcgisruntime.geometry.ProjectionConverter.convert - This seems to give me the option to suggest a wkid, x & y coord...  but It's poorly documented and, unfortunately, I can't get it to do any conversions)

    Dim wkid as integer= 102100
    Dim ref As New Esri.ArcGISRuntime.Geometry.SpatialReference(wkid)
    Dim converter = New Esri.ArcGISRuntime.Geometry.ProjectionConverter()
    Dim mappoint = New Esri.ArcGISRuntime.Geometry.MapPoint(-9638404.7886215933, 3990218.2857465576, ref)
    Dim typ As Type = GetType(Esri.ArcGISRuntime.Geometry.MapPoint)
    Dim glob As New Globalization.CultureInfo("en")
    Dim result = converter.Convert(mappoint, typ, wkid, glob)

FYI, if there are some arcgis account auth required, I do not have that.


texts: How to 'filter' a district by clicking on it in Leaflet?

texts:
I have created a basic layer of districts on a background google streetmap on QGIS. Using the qgis2leaflet plug in, I'm able to render a webpage of the same. I want the web product to be such that when a user clicks on a particular district (layer), it gets highlighted (to show that it is the district that has been 'selected') , zoomed in and displays the relevant information from the atribute table. How do I enable this filtering+zooming capability? (P.S: I read that the 'map events' funtionality in leaflet is the way to do it....but where exactly do I define these map events from QGIS desktop?)
Sorry if my question is stupid ......I'm new to geomapping softwares


texts: how to register view geometry made by MakeLine function?

texts:
I am trying to to create a spatial view on my spatialite database.

        CREATE VIEW v2_pump_view AS
          SELECT pump.ROWID AS ROWID, pump.id AS pump_id,
             MakeLine(
                 start_node.the_geom, end_node.the_geom
                 ) AS the_geom
          FROM v2_pump AS pump
          , v2_connection AS start_node
          , v2_connection AS end_node
          where pump.connection_start_id = start_node.id
          AND pump.connection_end_id = end_node.id
        ;

Then, to register with the views_geometry_columns:

    INSERT INTO views_geometry_columns
    (view_name, view_geometry, view_rowid, f_table_name, f_geometry_column)
    VALUES ('v2_pump_view', 'the_geom', 'ROWID',
        'v2_connection', 'the_geom');

Of course this doesn't work as the connection table has a point geometry. How I can register the line-geometry I made on the fly during the view creation by the MakeLine-function?


texts: How can I prevent KML files from being downloaded from Google My Maps?

texts: I have created a Map on Google My Maps using KML data from a third party company. I wish for anyone with the link to be able to view the map but I don't wish for them to have the option to download the KML data. Is there any way I can restrict downloading? Is there any encryption I can apply locally, or any other secure method?


texts: Creating cloud feature after applying Landsat8 Cloud Masking

texts:
There is a script in GEE to mask cloud:

// This example demonstrates the use of the pixel QA band to mask
// clouds in surface reflectance (SR) data.  It is suitable
// for use with any of the Landsat SR datasets.

// Load Landsat 8 surface reflectance data
var l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');

// Function to cloud mask from the Fmask band of Landsat 8 SR data.
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = ee.Number(2).pow(3).int();
  var cloudsBitMask = ee.Number(2).pow(5).int();

  // Get the pixel QA band.
  var qa = image.select('pixel_qa');

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));

  // Return the masked image, scaled to [0, 1].
  return image.updateMask(mask).divide(10000);
}

// Map the function over one year of data and take the median.
var composite = l8sr.filterDate('2016-01-01', '2016-12-31')
                    .map(maskL8sr)
                    .median();

// Display the results.
Map.addLayer(composite, {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.2});

After applying it, the cloud and shadow cloud pixels are labeled as "masked".
I want to create a feature class from the result of the script. In fact, I need to extract the masked pixels and make them as features.


texts: Spatial database containing coordinates for stars in the sky?

texts:
My goal is to map the Skydome ( stars in the sky projected in earth surface ) using GIS. Has anyone done this before?  I can't find any project close to this.

What I just need to know is if there is a dataset -not the paper or software, just the points (geometry) - of the relative position of the stars in the sky. I mean if I look at the sky I can see a kind of surface where the stars are. This is the Skydome. If this surface is over my table, I have a kind of "MAP" of the sky of my hemisphere. If I have both hemispheres I have the Sky Map of the stars.


texts:
I tried digitizing in to a new sqlite line layer with autoincrement id ticked.....

I want to create a 2-node line, but there seems to be a minimum limit of 3 nodes. Is this correct or am I doing something wrong? (left click to add a node, right click to complete the line)
once I right click the attribute dialog pops up no no auto incremented id.... I actually just want to digitize the lines with no further attribute data to be added and the id autoincremented. Have I missed something here?


texts: In terms of development, what capabilities does Google Earth Plugin have over Google Earth?


texts: Matching landbase data from different datasets

texts:
There is a layer containing all buildings from the electric landbase dataset.
And a layer containing buildings from the gas landbase.
The problem is, they are different and overlap each other. And I need to do the reconciliation of them both.

So, I do a selection of all the gas buildings which are WITHIN the 5m buffered electric buildings, and from that do a selection of all the electric buildings which are WITHIN the 5m buffered selected gas buildings.

For some building those shapes are very different between the electric and gas. So other method is to calculate the centroid points of all the buildings (gas and electric) and see how far from a gas building's centroid to the electric's. Together with comparison of shape area. If falling within an acceptable range then it's okay.

Arcgis 9 doesn't seem to have any tools supporting my second method.
Any other methods? Or how am I gonna do this?


texts: QGIS 3 Layers Panel - symbol color not the same as Categorized symbols on map

texts:
I have applied hex colors to the stroke and fill for a Simple Marker in the Symbol Selector from a join to an attribute table. Then I classified by a different column in the spatial table to give a text category in the Layers Panel. The colors on the map are correct, but the colors in the Layers Panel and the Layer Properties Symbology tab still show the color ramp. Is there any way around this?

example - the two points are the same classification - text removed


texts: Where does PostGIS store its topology tables?

texts: I have created a PostGIS topology table by importing a shapefile and following a tutorial online.  I can now load the edges and nodes into QGIS.  But I am not seeing the edges and nodes tables anywhere when I look at the database in PgAdmin.  Where are they being stored? I am making modifications to the edges and nodes tables in QGIS and I need to verify that they are being added properly.


texts: Coordinate System Assigmment and scale factors

texts:
I'm a bit unclear on what exactly assigning a coordinate system does in the model space of autocad (we are using 2015 civil 3d). Our local county coordinate systems have a scale factor associated to with the projection something along the lines of 1.00000534678....

Does that mean by assigning the coordinate system my drawing units change in the .dwg? As in 1 unit no longer equals 1 and now equals 1.00000534678....

I've done some tests and that does not appear to be the case. So, if not, what does assigning a coordinate system do? Why would it be beneficial for a land surveyor to assign a coordinate system in the .dwg? We pretty much have no reason to transform from one coordinate system to another. Other than using the automatic background map display (aerial imagery, road maps etc), I see no real reason to assign coordinates systems. Am I missing something?

We collect data in the field with GPS and start all projects with the county coordinate system assigned to the data collector so all of our coordinates are exported from collector into our .dwgs and are already in the proper coordinate system.


texts: Problem with encoding when clipping raster

texts:
I am using QGIS 2.18 on a Linux Mint computer using UTF8.
I downloaded solar radiation raster data from PVGIS. The files are encoded as us-ascii (according to the linux command tool file -bi).
Opening the data in QGIS works fine. But I want to clip this data to reduce it to the area of the project on which I work and here is the problem.

I created a .shp vector layer as the mask layer for clipping, and when I try clipping (Raster > Extraction > Clipper ) , I have the following Python error : Traceback (most recent call last):
File "/usr/share/qgis/python/plugins/GdalTools/tools/widgetPluginBase.py", line 88, in onFinished
self.base.onFinished(exitCode, status)
File "/usr/share/qgis/python/plugins/GdalTools/tools/dialogBase.py", line 208, in onFinished
msg = unicode(self.process.readAllStandardError())
UnicodeDecodeError: 'ascii' codec can't decode byte 0xc3 in position 52: ordinal not in range(128):

Obviously, it's something with the file encoding.
A quick search led me to this thread but that does not really help me solving the problem.


texts: Why Geoserver GWC REST put tiles into strange folders?

texts:
I make a GWC layer, and cache it using web interface. It puts the generated tiles in folders like:

gwc\temp3_layer_013D0001B7ED\EPSG_900913_11

After I want to reseed this layer using Geoserver REST Api. My request:

curl -v -u admin:geoserver -XPOST -H "Content-type: text/xml"
--data "<seedRequest>
          <name>temp3:layer_013D0001B7ED</name>
          <bounds>
            <coords>
              <double>6708811.26601302</double>
              <double>7677533.00995435</double>
              <double>6780820.92578184</double>
              <double>7753033.15577275</double>
            </coords>
          </bounds>
          <gridSetId>EPSG:900913</gridSetId>
          <zoomStart>10</zoomStart>
          <zoomStop>16</zoomStop>
          <format>image/png</format>
          <type>reseed</type>
          <threadCount>2</threadCount>
          <parameters>
            <entry>
              <string>STYLES</string>
              <string>style_013D0001B7ED</string>
            </entry>
          </parameters>
        </seedRequest>"
"http://localhost:8080/geoserver/gwc/rest/seed/temp3:layer_013D0001B7ED.xml"

But after reseed ends I do not see expected changes in gwc layer. In gwc folder I notice that REST request generate new folders like:

gwc\temp3_layer_013D0001B7ED\EPSG_900913_11_ae1eefd6f3b22d886816cc745224e69c74ebfde9

Why is this happening? And how to do reseed correctly?


texts:
I am working with the Natural Earth world map as a SHP file, but when I try to change the prime meridian (PRIMEM) of the coordinate system

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

I get this:

Projections seem to give the same results.  What is happening here?


texts: Simplification algorithms for 3D topology

texts:
My problem is about simplifying a given topology that can be opened ( an usual 3D reconstruction of a given square on a map ), totally closed (like a perfect sphere) or a topology with holes ( a sphere with a hole ).

Now, after reading a lot about this I probably should note that I already have a solution for simplifying the 3D mesh, my real problem is to keep the deviation and the difference at a really bare minimum level, in other words I have to preserve the topology in this process and at the same time lower the polycount.

The best that I found is the Douglas-Peucker algorithm, but it's for 2D shapes, I would like to hear a suggestion about an equivalent algorithm for 3D surfaces and topology.

Thanks.


texts: Splitting raster by polygons and defining output names using ArcMap

texts:
I need to split a raster file (DEM) by 2310 polygons (these are boundaries) using ArcGIS for Desktop. I used   the Split Raster tool but I have no control over the naming of the output rasters. They automatically get numbers assigned and I don't know which one is which. The shapefile table contains unique names and codes for each individual boundary (see below) and this is how I need to name my resulting rasters.
How can I split the raster and use a name/code from the polygon layer to name the outputs?
Either of them would work as long as I can later identify them. My ArcPy knowledge is very limited but if there's a solution that is not too complicated it could work too.


texts: Changing data and style with  in Leaflet popup

texts:
I'm loading a GeoJSON layer in Leaflet, and would like to be able to change the properties of each feature by using a  box in the feature's popup. When the feature's property changes, I'd like to change the style of the feature as well.

I'm sure this is possible- I've read several change style/property questions, and several javascript-in-popups questions, but just can't seem to get it right.

Here's a JS Fiddle: http://jsfiddle.net/quu3ykp6/7/

Desired workflow is as follows: user clicks a feature, changes the  box in the popup, this calls a javascript function that changes the style and properties of the feature in question.


texts: ENVI 5.0 Applying one ROI on various images

texts:
I'm comparing different classification methods (parallelepiped & max likehood) with ENVI 5.0 Classic. I'm supposed to use the same ROI from one image on various images and validate the classification with Kappa Index (Cohen's). How do I use the same ROI on different images? Note that I don't want to transfer the ROI from one image to another. What I need is to use the spectral information collected from ROIs on file A to train my classification methods for image B and C. The same training has to be used to different images. So far, I've tried using ROI Tool > calculate covariance with stats and then the Endmember Collection in order to apply the stats in a different image. I choose one of the methods described above and load my .sta file but the classification returns an entirely black image.

Any hints on how to solve this, or suggestions of a different method?

The images have been through atmospheric correction and cropping only.


texts: ArcGIS Online not saving data in a hosted feature layer

texts: I am trying to upload and save data in a hosted feature layer but every time I exit and reopen the layer my data is gone. I've tried saving it as a webmap and the data shows up on the webmap, but I need it to save in the feature layer (I am also trying to use this feature layer in a jupyter notebook). Any suggestions on what I should do?


texts: Calculate statistics for rasters with multiple bands (600) in QGIS

texts:
I am a new user of raster layers in QGIS.

I need to generate a new raster file with the average value for each pixel, from a raster file with 600 bands. The average has to be calculated from all the 600 bands, for each pixel.

The Raster Calculator only allows to add one band at a time, that would mean too much time and a risky process, considering that I have several similar raster layers that I have to work with.

Is there a way to add multiple bands at the same time in Raster Calculator in order to calculate it? Otherwise, is there another way to do it?

I tried also with RASS r.series but the output does not seem to be what I need, the layer is not shown in the canvas and the values are not what I expect to be (the data is about temperature)


texts:
I'm using Map Basic IDE Version 1.6.0.0. Every time I use the "." key the IDE always freezes. This freeze happens anytime I use this key, whether it be for comments, or to access a  variable's fields. Does anyone know how to fix this?
I'm assuming that using this key triggers some sort of auto-save or something but I am unsure.


texts:
Despite the topic Cannot open NetCDF File in QGIS, I cannot properly visualize my NetCDF file (nc4) in QGIS. It is quite strange because Panoply can do this (see image below). My NetCDF file visualize Standardized Precipitation Index (there are many results for the region of interest). It seems to me that I am able only to visualize it in some strange very narrow line if I loaded as raster.

This really doesn't help:

 Layer > Add Layer > Add Mesh Layer.

Error message:

Invalid Data Source: C:\Users\lukat\Downloads\IMERG_spi_gamma_03.nc is not a valid or recognized data source.


texts: Geographically weighted regression with QGIS / R / etc?

texts:
I have a small project where I'd like to investigate the relationship between number of events (many thousand x-y points) in a city and the rent of each neighbourhood (around 50 x-y centroid points of neighbourhood polygons, each with a corresponding rent value).

I have these two data sets in two csvs (one with all the events and their lat/long, the other with all the neighbourhood centroids and their rent value).

I was hoping I would be able to click some kind of magic button in QGIS to do a geographically weighted regression (I already did a basic non-GW regression and plot in R and there was a nice relationship) in order to account for the proximity of the actual points for events/neighbourhoods... but I can't see an option anywhere in the menu for GWR.

Therefore,

1) Is there a way to do this easily in QGIS? Perhaps a plugin of some sort?

2) If not, do I have to use R and what is the plugin for this?

I've never done GWR before.


texts: New expression of "legendInterface()" in qgis3.0?

texts:
What is the correct expression of "legendInterface()" in qgis3.0? I tried with "layerTreeNode", but it doesn't work.

root = QgsProject.instance().layerTreeRoot()

for group in [child for child in root.children() if child.nodeType() == 0]:
    if group=="POSTE_DETENTE":
        group_actif = root.findGroup(group)
        #now, I would iterate on layers in the list L in order to get the polygon layers
        for (type, sys, champ, name) in L:
          #get the polygon layer
          if type=='polygon':
            #create_layer is a function that I create previously and return a new layer in memory
            layer=create_layer(type, sys, champ, name):
            # add layer to the legend
            QgsProject.instance().addMapLayer(layer, False)
            # add layer to the group - False to let the choice where to place the layer in the group
            group_actif.addLayer(layer)


texts: How to calculate the intersected area from a layer with overlapping buffers in QGIS

texts:
I need to calculate the area of landusetypes in 50m around various dots (all in one shapefile). My problem is, the buffers around the dots are overlapping each over in some places. I tried to intersect/clip the layer (e.g. grassland) with the layer "buffer", but the overlapping parts are only counted once to the overlapping buffers not twice.

Is there any tool in QGIS which might help me? Or is there any other way to get this output?


texts: Google Map Projection Skew - Need correct Projection

texts: I'm having a really difficult time resolving skewed grids when projecting onto Google maps using Alteryx outputs. Within the module, I'm taking GRD files and aggregating rows and columns to re-size from 30m to 250m grids. Then projecting using EPSG:3857 which I understand to be the correct Web Mercator Aux Sphere for Google/Bing. However, no matter what I research and attempt to implement, I always end up with skew in my squares. Any suggestions would be very helpful. Thanks


texts: Get all locations of Panoramio photos in a given extent

texts:
To create a visual quality heat map I'm trying to get a list of all Panoramio photos in an about 2000 Km2 wide study area. My idea is to perform a density analysis in GIS on a point shapefile including all the photo locations. I've used the well-known Panoramio syntax to get the coordinates:

http://www.panoramio.com/map/get_panoramas.php?set=full&from=0&to=499&minx=12.068&miny=42.921&maxx=12.577&maxy=43.371&size=medium&mapfilter=false

But I've encountered two kind of problems:

Every single output seems limited to 498 records (not only 100 as declared by Panoramio..)
The number of photos seems reduced automatically depending on the bounding rectangle extents, so not all the photos within the area are returned even setting consecutive from-to ranges.

So a possible solution could be to loop the query according to a moving window approach (let's say 0.01x0.01 deg wide) and to merge the results.

Unfortunately I don't have experience on programming...
Any suggestion about how to do this?


texts: Translate Orientation from Sensor data to Geo location

texts:
I am tracking a device that might loose GPS signal then and there. And when the GPS signal is too low, I should use the accelerometer and Gyroscope sensors data to guess which way they should have gone.

I am not doing it in real time. I can store the data and send up for post processing.
So all I want to do is from the point (Lat1, Long1) where we had  a good GPS, I need to find new location (Lat2, Long2) from the sensors data.

I was looking at different algorithms. One of them was inertial Navigational System. I applied Kalman filter and I was able to get to a point to find the Quaternion for each data.

Now can some one help me to translate the orientation data to an actual position? What algorithm or formula should I use? Or is there any program available in C# or in any other language?


texts: Reclass Using Field Calculator in ArcMap

texts:
I'm trying to reclass 3000, 2000, 1000 and 0.5 into 1, 2, 3 and 4 respectively in a new field using field calculator. However it keeps returning an error message and I can't figure out why, any ideas?

def Reclass( !BUFF_DIST! ):
    if !BUFF_DIST! = 3000:
        return 1
    else if (!BUFF_DIST!) = 2000:
        return 2
    else if (!BUFF_DIST!) = 2000:
        return 3
    else:
        return 4


texts: Is there anything obviously wrong with this hit test code?

texts:
I want to know if the IPoint of the mouse position is hitting an IPolyline, I wrote this code, and was optimistic that it would work, but it always returns false.

  IPolyline poly1 = pc1 as IPolyline;

        IHitTest test = poly1 as IHitTest;
        if (test != null)
        {
            int a = 0;
            int b = 0;
            bool c = false;
            double distance = 0;
            IPoint hitPoint = new PointClass();
            bool isHit = test.HitTest(mousepoint,
                                    50000,
                                    esriGeometryHitPartType.esriGeometryPartNone,
                                    hitPoint,
                                    ref distance,
                                    ref a, ref b, ref c);

        }


texts: Subsetting single raster from raster list using ArcPy?

texts:
I have a list of rasters, created by arcpy using ListRasters. Rasters are stored in a geodatabase. How can I subset only one raster of the raster list? The goal is to make a new list, but containing only one raster file.

Ultimate goal is to loop through list of rasters and calculate Area solar radiation for each one of them. But, but for learning/training purposes I wish to explore my parameters only on single raster for now. (I know I can read-in just single raster, but I thought that subsetting might be easy as well?)

I have thought that I can simply subset the element by the index value, but this seems not working in arcpy. Any advices? I think I am missing something very obvious...

import arcpy, os

# Set working environment
arcpy.env.workspace = os.path.join(inWD, "output/bufRastTwins.gdb")

# List all rasters with defined name wildcard: r_1, r_2, r_3.. to r_100
allRasters = arcpy.ListRasters("r_*")

# Select the first raster from the rasters list
myRasters = allRasters[0]
print myRasters

My attempt leads to loop through raster name, not throught the subsetted raster itself:

print myRaster:
r_1  # first element succesfully created

but for raster in myRasters:
   print raster    # this is obviously wrong

r
_
1


texts: Help with field calculator in Python script

texts:
I am getting a syntax error in my script but I know from testing it in ArcMap attribute table, the syntax is fine. Anyone have any idea why this code doesn't work? I see that a search cursor has been used in this situation but (a) I'm not familiar with them and (b) why wont the way I'm doing it work?

EDIT* - I have added an alternate working codeblock to the script below. Seems both ways work just fine!

import arcpy, os
from arcpy import env

#Path of the workspace
in_workspace = r"C:\folder"

fieldName = "Symbology"

expression = "Class(!City!)"

This is also another way that it works. So both work just fine. Thanks again for all the help!

code ='''
def Class(City):
    if 'Milwaukee' in City:
        return 2
    elif 'Kenosha' in City:
        return 4
    elif 'Madison' in City:
        return 9
    else:
        return 0
'''

for dirpath, dirnames, filenames in arcpy.da.Walk(in_workspace, datatype="FeatureClass",type="All"):
    if "skipFolder" in dirnames:
        dirnames.remove('skipFolder')
    for filename in filenames:
        if "City" in [f.name for f in arcpy.ListFields(os.path.join(dirpath, filename))]:
            arcpy.AddField_management(os.path.join(dirpath,filename), fieldName, "SHORT", "", "", "", "", "NULLABLE")
            arcpy.CalculateField_management(os.path.join(dirpath,filename), "Symbology", expression, "PYTHON", code)
        else:
            print "no City field for dataset ", filename

#Prints text that shows script has finished running
print "message that reflected field calc use"
print "Finished script"


texts: How to display only one pie chart for polygon features in QGIS?

texts:
I'm working with Qgis on a project and I have a small problem.
I have to do some demographic research and therefore I have a .shp file from a given country with 10 districts for example.

Now, I know how to insert pie charts but the problem is that I only want to insert one pie chart for one district because this area is of importance to me for example.

So I added a column in the attribute table with data for the district that I want to display and left all the other values 'NULL'.

However, if I insert pie charts, they appear for all the districts, even though I've not inserted any values in the attribute table for the other 9 districts that are of no interest to me.

Does anybody know how I can only display one pie chart for 1 feature of the polygon instead of having 1 good pie chart and 9 empty on the map?


texts: Importing classification bitmap made by another software into QGIS?

texts:
I am currently trying to import a classification bitmap made by our own deep learning software into QGIS. The classification is saved in a TIF file, the classes are 8-bit integers. We have the class names in a text file and the coordinates (MTM) to geo-reference the image in another text file.

Ideally, we would want QGIS to assign colours to the classes, and when the user clicks on a particular pixel of the image, to have QGIS display the class name.


texts:
Here is my code so far:

import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

workspace = arcpy.env.workspace = "C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6"
arcpy.CreateFolder_management(workspace, "NDVI")
arcpy.env.overwriteOutput = True

RED_1 = arcpy.Raster("C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6/landtm1.img/Layer_3")
NIR_1 = arcpy.Raster("C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6/landtm1.img/Layer_4")

RED_2 = arcpy.Raster("C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6/landtm2.img/Layer_3")
NIR_2 = arcpy.Raster("C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6/landtm2.img/Layer_4")

RED_3 = arcpy.Raster("C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6/landtm3.img/Layer_3")
NIR_3 = arcpy.Raster("C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6/landtm3.img/Layer_4")

print "Program recognizes raster bands."

NDVI_1 = (NIR_1 - RED_1)/(NIR_1 + RED_1)

sr = arcpy.sa.Con(RED_1<>0, NDVI_1, -9999)

NDVI_2 = (NIR_2 - RED_2)/(NIR_2 + RED_2)

sr = arcpy.sa.Con(RED_2<>0, NDVI_2, -9999)

NDVI_3 = (NIR_3 - RED_3)/(NIR_3 + RED_3)

sr = arcpy.sa.Con(RED_3<>0, NDVI_3, -9999)

print "NDVIs passed"

arcpy.CopyRaster_management("NDVI_1", "C:/Users/KJanus/Desktop/Intro_ArcGIS/RasterData_Lab6/NDVI",
format= "TIFF")

When I run the code I get these two errors: ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Raster: Dataset NDVI_1 does not exist or is not supported


texts:
Does MongoDB follow OGC compliancy standards for GeoJSON objects, or some other standard?

I have not been able to find a direct statement on this. A follow up to this is whether there are any tools/methods for taking a GeoJSON object and converting it to be compliant with the standard that Mongo uses?


texts: How to define Z dimension on table with PostGIS for import in TileMill?

texts:
I'm trying to reproduce the idea behind this blogpost from Mapbox, but with other data, for example the ne_10m_populated_places_simple dataset from Natural Earth.  By adding a z-dimension to the data (eg. the scalerank column, size of city), the idea is to be able to set a height in CartoCSS in Tilemill and get a 3D-representation of the data, with x,y as geo-point and z as value (like in Mapbox's example).

The result could perhaps be something like this (except this map is oblique, different story):

What do I need to do with PostGIS (or other) so that the z-dimension is set for value and correctly converted so I can set the height in CartoCSS?


texts: Converting .shp file to .csv file using ogr2ogr?

texts:
I want to covert the .shp file to .csv file by using ogr2ogr
and I am running

PGCLIENTENCODING=UTF8 ogr2ogr -f CSV "/home/pj/Desktop/test1.csv" "/home/pj/Desktop/20160705102900718.shp"

but the shell prompt

ERROR 4: Unable to open /home/pj/Desktop/20160705102900718.shx or /home/pj/Desktop/20160705102900718.SHX.
FAILURE:
Unable to open datasource `/home/pj/Desktop/20160705102900718.shp' with the following drivers.
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> SDTS
  -> TIGER
  -> S57
  -> DGN
  -> VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> NAS
  -> GML
  -> GPX
  -> LIBKML
  -> KML
  -> GeoJSON
  -> Interlis 1
  -> Interlis 2
  -> GMT
  -> SQLite
  -> DODS
  -> ODBC
  -> PGeo
  -> MSSQLSpatial
  -> OGDI
  -> PostgreSQL
  -> MySQL
  -> PCIDSK
  -> XPlane
  -> AVCBin
  -> AVCE00
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDump
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> PDS
  -> WFS
  -> HTF
  -> AeronavFAA
  -> Geomedia
  -> EDIGEO
  -> GFT
  -> SVG
  -> CouchDB
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> XLS
  -> ODS
  -> XLSX
  -> ElasticSearch
  -> PDF
]


texts: How do I create a mapserver config file for tiles downloaded from WMTS service?

texts:
I have downloaded Austria official map tiles from http://www.basemap.at/index_en.html using the provided offline downloader tool. Tiles are stored in my hard disk with paths similar to /<level>/<number>/<number>.png .

I now need to create a .map configuration file for mapserver, however all I have is the GetCapabilities result from the provider: http://www.basemap.at/wmts/1.0.0/WMTSCapabilities.xml .

Being new to GIS, I am unsure about whether I can manually create a MapServer .map file from scratch based on the GetCapabilities information, or if there is any tool or approach I shall use.

What is the preferred approach to create a .map file for a given set of downloaded raster tiles?


texts: Connecting line features and determining the length of the longest line

texts:
I have a line feature (see image) representing a river that I created using the Stream_to_Feature tool. The attribute table contains several records representing different lines - the problem is the longest line (easily distinguishable visually) is not represented as a single line in the table, it is actually made up of many smaller lines. The lines appear to be touching, although they do not cross through each other.

How can I merge these lines and then determine the length of the longest line using ArcObjects or manual methods that I can convert to ArcObjects? An even better solution would involve getting rid of all of the tributaries and leaving me with just the river channel as a line.


texts:
I am trying to import an OSM file of a State Park that was added in OpenStreetMap into QGIS, and it seems no matter what I do, I just cannot get it to work.

I've installed a bunch of plugins, etc, but nothing.


texts: how to copy layer and add to same mxd and rename layer using ArcPy and ArcMap 10.2?

texts:
I am trying to copy a layer in an mxd and add that layer to the bottom of the TOC using Python, but I have not found a way to save the mxd programatically without making a copy.  I also haven't found how to rename the copied layer.
Here is my code: it runs with no errors but the copied layer is not added to the mxd when I open it.

import arcpy
import arcpy.mapping

mxd = arcpy.mapping.MapDocument(r"D:/arcgisserver/directories/my_mxd.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
print "df : " + df.name
addLayer = arcpy.mapping.ListLayers(mxd, "", df)[0]
print "layer ? : " + addLayer.name

arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")

# Refresh things
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd, df, addLayer


texts: Identifying counties further than ten miles from bank using ArcGIS Desktop?

texts:
If I have the latitude and longitude of all US banks.

How can I identify counties whose center is ten miles or further from any bank using that data?


texts: Accessing Vector Field Parameter Value inside Python Script

texts:
I have created a QGIS Graphic Model using the built-in Graphic Modeler. My model is supposed to ask the user for an input layer, an overlay layer, and a field from the overlay layer.

The model then should calculate the area-weighted average of that field for each feature in the input layer. The result will be an input layer with an extra field for the area-weighted average. I have the complete workflow in my mind. Till now I have reached this point, where I want to multiply the field given by the user with the area of the feature in the intersection layer.

I have converted the model to a python script to use the value of the parameter 'Field to Average' in naming my field in the 'Field Calculator' and in calculating the new field. My question is how do I access the value of Field to Average parameter inside my Field Calculator Formula and Field Name. I have tried using parameter(fieldtoaverge) and it is not working. See my comments in the script inside Field Calculator Algorithm

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsProcessingParameterField
import processing

class Model(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('overlaylayer', 'Overlay Layer', types=[QgsProcessing.TypeVectorPolygon], defaultValue=None))
        self.addParameter(QgsProcessingParameterVectorLayer('inputlayer', 'Input Layer', types=[QgsProcessing.TypeVectorPolygon], defaultValue=None))
        self.addParameter(QgsProcessingParameterField('fieldtoaverage', 'Field to Average', type=QgsProcessingParameterField.Numeric, parentLayerParameterName='overlaylayer', allowMultiple=False, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(2, model_feedback)
        results = {}
        outputs = {}

        # Intersection
        alg_params = {
            'INPUT': parameters['inputlayer'],
            'INPUT_FIELDS': None,
            'OVERLAY': parameters['inputlayer'],
            'OVERLAY_FIELDS': None,
            'OVERLAY_FIELDS_PREFIX': '',
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['Intersection'] = processing.run('native:intersection', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(1)
        if feedback.isCanceled():
            return {}

        # Field calculator
        alg_params = {
            'FIELD_LENGTH': 10,
            'FIELD_NAME': 'Field Avg * Area', # I need to name it with 'field selected by user * Area'
            'FIELD_PRECISION': 3,
            'FIELD_TYPE': 0,
            'FORMULA': ' parameter(fieldtoaverage) *  $area ', # This is not working, how do I acess the field selected by the user here
            'INPUT': outputs['Intersection']['OUTPUT'],
            'NEW_FIELD': True,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['FieldCalculator'] = processing.run('qgis:fieldcalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        return results

    def name(self):
        return 'model'

    def displayName(self):
        return 'model'

    def group(self):
        return ''

    def groupId(self):
        return ''

    def createInstance(self):
        return Model()


texts:
I am trying to create WMS service using Postgres data in MapServer, my map file (wms.map) is:

MAP
  EXTENT 8240841.450745 2400319.900386 9218007.463612 3107183.406947
  TRANSPARENT on
  IMAGETYPE png
  size 150 500
  web
    metadata
      "wms_title" "Customer Airports"
      "wms_onlineresource" "http://localhost:8085/cgi-bin/mapserv.exe?"
      "wms_enable_request" "*"
      "wms_srs" "EPSG:4326 EPSG:900913"
      "wms_feature_info_mime_type" "text/html"
      "wms_format" "image/png"
    end
  end
  LAYER
      CONNECTIONTYPE postgis
      CONNECTION "user=postgres password=abcd dbname=sample port=5432 host=localhost"
      DATA "geom FROM district_old using unique gid"
      STATUS Default
      TYPE polygon
      PROJECTION
         "init=epsg:900913"
      END
      labelitem 'district'
      class
          style
              outlinecolor 90 90 90
              width 2
          end
          LABEL
             size 2
             BACKGROUNDCOLOR 255 201 14
             position auto
          END
      end
   END
END # nmtoken edit added missing Map END tag

When I paste in the browser this URL:

"http://localhost/cgi-bin/mapserv.exe?map=c://ms4w/apps/webgis/wms.map&SERVICE=WMS&REQUEST=GetCapabilities" it is not displaying the XML file instead it is downloading a mapserv.exe file

Please correct my code.


texts:
I am trying to create a fishnet of x by x miles for all USA, same as this basically:
http://www.giscourse.com/creating-a-fishnet-grid-using-arcgis-10/

However, I want to do this with python.

My current code works but it doesn't take distance in miles, but rather in degrees, and therefore it doesn't take into consideration how for different lat, lon points the conversion of distance to degrees changes (Calculating longitude length in miles?).

This is my current code:

from pygeoif import geometry
import shapefile
import numpy as np

sq_len_deg = 1
lats = np.arange(26, 50.1, sq_len_deg)
lons = np.arange(-125, -69.9, sq_len_deg)
a=[]
geometries=[]
for lat in lats:
    for lon in lons:
        a.append([(lon,lat),(lon,lat-sq_len_deg),(lon-sq_len_deg,lat-sq_len_deg),(lon-sq_len_deg,lat)])
        geometries.append('POLYGON ((%s %s, %s %s, %s %s, %s %s))'%(lon,lat,lon,lat-sq_len_deg,lon-sq_len_deg,lat-sq_len_deg,lon-sq_len_deg,lat))
w = shapefile.Writer(shapefile.POLYGON)
w.field('test','N')
for i,poly in enumerate(geometries):
    w.poly([geometry.from_wkt(poly).exterior.coords])
    print(i)
    w.record(i)
w.save('officialpysh')

I can see this going into two directions, some function already exist and I can use this, or a suggestion on how to modify this code.


texts: spatialite - azimuth function not available

texts:
I am working with last spatialite versions on both Windows (7) and Ubuntu (13.10).

I installed these from reference spatialite website. For both, version that I downloaded is 4.1.1 (nevertheless in Ubuntu, at spatialite start up, i get a version number : 4.1.0 ; I assume this is a little bug)

Details on the packages :

binaries for windows: spatialite-4.1.1-win-x86.zip
source for ubuntu: libspatialite-4.1.1.tar.gz and all other usefull ressources

I then developped a sql script on windows, notably using azimuth /st_azimuth
this works fine on this OS. I have problems when I use it on Ubuntu and get the following error:

no such function: Azimuth

after some search, I undesrtood that this could be related to a missing lib: liblwgeom.

I had some more search to get this lib installed properly (at least no error message at the end) using the procedure described here (I effectively did the make and make install), but this does not seem to help.

Questions:

is the "wrong" version number for spatialite 4.1.1 a known bug?
is the missing Azimuth a know problem (i could not find any info)
is installing liblwgeom the solution to my probleùm ?
how to know if the liblwgeom install properly working ?
what may I have missed in the process ?

Thanks for your help and pieces of advice


texts: Mapping GPS coordinate from Google Maps in QGIS

texts:
I have a number of addresses that I would like to map as point data (along with attributes) in QGIS. I got the GPS coordinates from Google Maps, but when I create the layer from the table, the points are way off. I have just copied and pasted the GPS coordinates from Google maps into a spreadsheet. See image below.

They should be located around Salinas, CA, USA not off the coast of Africa. I have tried using different Coordinates Reference Systems, for example: Google Mercator, NAD83(CSRS), NAD83, WGS 84 / Simple Mercator, WGS 84 / Pseudo Mercator. I am using the Google Maps base layer from the open layer plugin, yet none of these project the points where I would like them to be.

Is it the format of the GPS coordinates?

Or I am totally using the wrong Coordinates Reference System?


texts:
Running on Mac OS.

I am using a composer template that I have used in the past with this version of QGIS (3.2) that uses an atlas. Now it does not display the Atlas tab in the composer.  The generate atlas box is checked.

I then installed 3.4 with the same result.
Next I tried opening an empty composer and loaded a map, still no Atlas tab and checking the generate atlas box had no effect (it did not register -- no tick -- it appears disabled).

I swapped to a near identical system and started QGIS 3.2 opening the same project (it is stored on Google Drive) using the same PostGIS database with the same result.

Anyone have any idea what is going on?   I last used this print template a few months ago.


texts: Creating different label styles for different polygon classes using QGIS?

texts: I'm new to QGIS 2. I have a polygon shape with two classes, and I want to set different labeling style and different "Scale-based visibility" but I don't get how the expression works.


texts: Removing black background and black font from raster in QGIS

texts:
I have a raster where I want to remove the black background.

I tried it with setting "Additional no data value" in Layer properties -> Transparency to 0.

But the place names on the raster are also black and they disappear as well. Is there another method that would work here?


texts: Generate raster from factor-linked data in R

texts:
I have a land classification raster with factor levels set (pixel val of 1 = 'type a', etc), plus a data frame holding the same factor in the first column, and then several columns of numeric attribute data.

I'd like to generate a new rasterBrick/rasterStack where that numeric data is mapped out - one map for each attribute. Not sure how to do this properly in R. I've tried

test <- merge(raster@data@attributes[[1]], data.frame, by = 'category')

but that gives me a single-band rasterLayer with multiple attributes (i.e. a RAT with more columns), which I can't then write to GTiff properly.


texts: How to add pixels offset based on resolution

texts:
Currently I'm displaying popups on my map using this utility by Matt Walker.

My code is something like this:

map.defaultPopup = new ol.Overlay.Popup({ panMapIfOutOfView: false });
map.addOverlay(map.defaultPopup);

// ...

var select = new ol.interaction.Select(options));

select.getFeatures().on('change:length', function (e) {
    this.map.defaultPopup.hide();

    if (e.target.getArray().length > 0) {
        var feature = e.target.item(0);
        var center = feature.getCenter(); // getCenter is defined elsewhere
        //center[1 /* Y */] += ????;
        map.defaultPopup.show(center, "
Lorem ipsum dolor sit amet.
");
    }
});

map.addInteraction(select);

I can add an offset "????" and it works, but it is not resolution independent (the popup gets too far away from the feature on higher zoom levels).

I want to add a value that moves the location of the popup by the desired amount of pixels, regardless of the resolution. Is there any way to achieve this?


texts:
I am trying to leverage my GeoServer and add a WMS layer to a map created with of Google Maps.

What I am getting instead is the tile repeated across the display. I am passing the BBox from the current map extent to the WMS request which works properly as a standalone URL but somehow displays incorrectly on the map.

I am familiar with GeoServer and the WMS request but new to Google Maps JavaScript API.

Here is what I have tried and the resulting image.

    google.maps.event.addListener(map, "bounds_changed", function () {
        var censusLayer =
            new google.maps.ImageMapType(
                {

                    getTileUrl: function (coord, zoom) {
                        var ne = map.getBounds().getNorthEast();
                        var sw = map.getBounds().getSouthWest();
                        var bbox = sw.lng() + "," + sw.lat() + "," + ne.lng() + "," + ne.lat();

                        var container = document.getElementById('map');
                        // from geoserver
                        //http://104.215.74.199:8080/geoserver/hcad/wms?service=WMS&version=1.1.0&request=GetMap&layers=hcad:parcels_simple_4326&styles=&bbox=-95.96082688782865,29.497301702948995,-94.90761061770445,30.170549317160965&width=1538&height=517&srs=EPSG:4326&format=image%2Fpng
                        // from application
                        //http://104.215.74.199:8080/geoserver/hcad/wms?service=WMS&version=1.1.0&request=GetMap&layers=hcad:parcels_simple_4326&styles=&bbox=-97.44697158203127,29.149074739047958,-93.22272841796877,30.381673782887290&width=1538&height517&srs=EPSG:4326&format=image%2Fpng

                        //base WMS URL
                        var url = "http://104.215.74.199:8080/geoserver/hcad/wms?"

                        url += "service=WMS";           //WMS service
                        url += "&version=1.1.0";         //WMS version
                        url += "&request=GetMap";        //WMS operation
                        url += "&layers=hcad:parcels_simple_4326"; //WMS layers to draw
                        url += "&styles=";               //use default style
                        url += "&bbox=" + bbox;          //set bounding box for tile
                        url += "&width=256" //+ container.offsetWidth;
                        url += "&height=256" //+ container.offsetHeight;
                        url += "&srs=EPSG:4326";         //projection WGS84;
                        url += "&format=image%2Fpng";      //image format
                        url += "&TRANSPARENT=TRUE";      //only draw areas where we have data
                        url += "&tiled=true";

                        console.log(url);
                        return url; //return WMS URL for the tile
                        return container;
                    },

                    tileSize: new google.maps.Size(256, 256),
                    opacity: 0.5,
                    isPng: true
                });

        //add WMS layer to map
        map.overlayMapTypes.push(censusLayer);
    });
}

And the result so far...


texts: How to load/enable r.stream.distance and r.lfp in QGIS 3.6.1 with GRASS 7.6.1 (Windows 10)?

texts:
I am new to QGIS and GRASS. I would like to use r.lfp in QGIS with GRASS.

So far I installed addons extensions r.stream.distance and r.lfp in GRASS (Settings/Addons extension/Install extension from addons [g.extension]), without checking the two optional check boxes (system-wide and verbose).

I was hoping opening QGIS with GRASS after that would allow me to see the two new raster commands in QGIS with GRASS (Processing Toolbox/GRASS/Raster (r.*)). But this is not the case.

Any idea what I am missing and how to correct it?


texts: Resolution of results exported from Google Earth Engine

texts:
I use this code to export my classification result from Sentinel2 images in GEE:
// Export the results to Drive
Export.image.toDrive({
  image: result,
  description: 'Unsupervised Classified Image',
  scale: 20,
  crs: EPSG=32600,
  region: geometry
});

And I also choose 10m as scale while exporting to Google Drive.
I thought I am selecting the output resolution by setting these vlaues to 10.
When I check the information of the exported tiff in QGIS, the Pixel Size is:   [8.983152841195214829e-05,-8.983152841195214829e-05].
How can I obtain the 10m resolution in my output?


texts: measuring bankfull width by distance between points

texts:
I have a large set of polygons that represent sections of streams, each of which has a Unique ID. Each of these polygons is made up of a series of left bank points (LB1,LB2,etc.) and right bank points (RB1,RB2,etc.). I need to calculate the average bankfull width (width of polygon) for these streams.

I think something like point distance in arcmap would work, but instead of finding the distance based on nearest point I want to match up LB1 and measure the distance to RB1 and sequentially do this for each set up bank points.


texts: How to draw a rectangle graphic element (given feature and extents in meters) with ArcObjects?

texts:
I am trying to render a rectangle graphics box onto the map, given the center location of an polygonal feature and the extent dimensions (length and width) in meters. The feature has the spatial reference of its shape set to WGS84 though.

I am able to get the location of the feature (an IPoint) by taking its centroid. Is there then a straightforward way to accomplish the drawing of such a box with meter-dimensions around that point?

Looking at the API, IEnvelope seems to be what I need, but will I have issues with unit-conversions (since I need the extent dimensions to be in meters)? And how should I go about rendering the IEnvelope?

Any help is appreciated, thanks!


texts: Mapinfo Vertical Mapper interpolation VS QGIS

texts: I am currently trying to fade out the use of MapInfo and only use QGIS (or Global Mapper) however I find that I am still needing to use Vertical Mapper for interpolations. We are using MapInfo V12 and QGIS V2.8. Below are images of Vertical Mapper grid interpolation using Nearest Neighbour and QGIS grid interpolation using nearest neighbour. Vertical Mapper's output is much more smooth and allows for smooth contour lines to be created. I guess what I'm asking is can QGIS produce these results?


texts: Why would a new install (ArcGIS Desktop 10.4.1) not be able to open existing MXDs without crashing?

texts:
A coworker can open up a brand new install of 10.4.1, can add data, and save a map, but cannot open existing MXDs. The system tries for a while, then shuts ArcMap down with no error.

Is there any chance this has to do with templates?

Currently, due to some profile issues, I have limited access to his machine for testing.

OTW, it might take days for an administrator to do a new install. I verified that Windows is not hiding files or file types in the template directory on his C drive.


texts:
I have download the following Landsat 7 image

LE07_L1TP_231095_20011014_20170202_01_T1.

I want to open a true color RGB composition using bands 1, 2 and 3. For that, I import the Tiff files of the bands. When a create the RGB composition (raster, miscellaneous, virtual raster) the result does not look very nice. My issue is not that the image is dark (there is a lot of ice/snow) but that the true colors are not very accurate. Water looks purple, glaciers have some purple areas, etc.
Am I missing something to create a good looking true color RGB?

The output should look like this

-- EDIT --


texts: How to create routable network by pgrouting with assign_vertex_id?

texts: I'm trying to create the routable network as per http://underdark.wordpress.com/2011/02/07/a-beginners-guide-to-pgrouting/and use assign_vertex_id. Do I still need to create the node table? How would code for creating network look like?


texts: Excluding points based on several features using ArcGIS ModelBuilder?

texts:
I have 3 separate feature classes:

Lines: 3 pcs, do not touch eachother

Points: 6pcs, at the start and end vertex of the lines

Polygons: 1 pcs, derived from joined buffers around the points based on the top line.

I only have the geodata, no additional fields of ID's to connect the features to eachother. I want to create a new feature class with just the two points that are both within the polygon and are based on the top line.

The 2 top lines are separate even though they look linked in this resolution.

I would like the solution to be able to be reproduced in the model builder as I need to iterate this method over many scenarios.

The polygon is one feature, not 2 as I originally stated. They are a union from two previously separated buffers


texts: QGIS 3: How to count the mean of shapefile attributes that lie in a specific area

texts:
so I have a shapefile with buildings and the shapefile has the height of the buildings in the attribute table. I have a second shapefile with the districts the buildings are located in. Now I would like to count the mean height of the buildings in each district and add it to the attribute table of the districts.
For raster data it would be easy using zonal statistics to get the mean for a specific area, but I cant figure out the solution for 'shapes within shapes'.


texts: Finding features with equal distance for houses in archaeological data using QGIS

texts:
My goal is to find house sketches in archaeological data.
I have a shapefile from an excavation where all the archaeological features, mainly pits are measured as polygons. The polygons have different sizes and shapes, but are mainly roundish.
The assumption is, that houses in prehistory consist of pits that stand in line and have roughly the same distance to another or make up a rectangle or are symetrically arranged.
Is there a way to highlight pits/features in the data that match the above criteria so that e.g.

All the features that build up a rectangle or construct a right angle are marked or the rectangles/right angles are drawn to another layer
All the features that are symetrically arranged are selected or marked in the same colour
All the features where the distances inbetween the features are similar (e.g. in between a range form 180 cm to 220 cm) are selected?

EDIT 23.04.21
So i cleaned the Data and selectet some features that are likely to be Postholes. After that i started to change the features with the tool centroid to points.
The idea is now to connect all the dots with lines and select the Houses via the lenth of the lines. I tried to connect the dots buy using the delaunay triangulation and convert the polygons into lines afterwards. Sadly i dont get only single lines but also i have lines that have still a triangular form. Is it possible to convert the polygons in a way that i have only lines without the vertices?
Thx in advance!


texts: QGIS Spatialite tabular updates are very slow

texts:
I am exploring Spatialite in QGIS (2.10.1), and it is very slow to save tabular edits (with 1200 records, several minutes once I click the pencil button, versus less than 1 second with the data as a shapefile!!!).  Here's my workflow:

Add Vector a shapefile into QGIS and then Save As a Spatialite database.
Add Vector the Spatialite layer from above and Open Attribute Table.
Toggle Edit the attribute table and add an integer field, calculating it.
Toggle editing off and save.
File saving takes 2-3 minutes while the hard drive is on continuously.

When I repeat the above workflow using the same data in shapefile format, the file saving takes less than 1 second.

What can I do to speed up Spatialite?


texts: gdal_translate slightly changes output pixel values

texts:
I have used this batch script

for /r %%g in (*.jp2) do gdal_translate -co COMPRESS=LZW -co
NUM_THREADS=ALL_CPUS %%g %%~dpng.tif

to recursively convert all my .jp2 files in .tif files. However, the output .tif file doesn't have exact pixel values as the original .jp2 image (see in the picture).

On the left side is an example on a 10 meter resolution raster, and on the right the same problem is present in the 20 meter raster.

How can i correct this issue?

I need the tif files to have the exact same values as the original.


texts: pgRouting only partially segments network?

texts:
I have a road network traced from an old map with 360 MultiLineString features in table aroads. Many meet at endpoints or intermediate vertices, but some intersect.

When I run SELECT pgr_nodeNetwork('aroads', 0.0001) many are segmented properly (now 899 features in aroads_noded), but many are skipped. pgr_analyzegraph() detects 79 intersections before and 27 after.

Can't split this intersection (and others) for the life of me, any ideas?

Before (aroads) --

NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 33
NOTICE:                            Dead ends: 193
NOTICE:  Potential gaps found near dead ends: 116
NOTICE:               Intersections detected: 79
NOTICE:                      Ring geometries: 0

After (aroads_noded) --

NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 6
NOTICE:                            Dead ends: 78
NOTICE:  Potential gaps found near dead ends: 1
NOTICE:               Intersections detected: 27
NOTICE:                      Ring geometries: 0


texts:
I've unfortunately experienced a version control issue between offices where one person was editing an older version of a polygon shapefile. Now I need to reconcile those edits into the current version. I've noticed that ESRI has a great tool for identifying polyline changes between shapefiles with the Detect Features Changes tool. However, it doesn't appear that a similar tool exists for polygons.

Any suggests for how I can quickly identify which polygons have been edited? Keep in mind, there are over a thousand polygons. And no, I am not on an enterprise GIS with versioning.