texts: Error Loading Layer/Shapefile in a Standalone Python QGIS Application

texts:
I tried loading a shapefile on Python using PyQgis API but to no avail. I double-checked the path to shapefile and found it to be correct. QGIS module also seems to be imported just fine. When I checked on the provider list in QgsRegistry, it returns nothing.

May I know what I'm missing or how I should troubleshoot?

I am using Ubuntu 12.04, QGIS 2.4.0 Chugiak and Python 2.7.3.

Following are my output and code:

  /usr/bin/python2.7
  /home/victorzhiyulee/IdeaProjects/Delineation/select_dun_calculate_print.py
  Application state: QGIS_PREFIX_PATH env var: Prefix: /usr/bin/qgis
  Plugin Path: /usr/bin/qgis/lib/qgis/plugins Package Data Path:
  /usr/bin/qgis/share/qgis Active Theme Name: Active Theme Path:
  :/images/themes// Default Theme Path: :/images/themes/default/ SVG
  Search Paths: /usr/bin/qgis/share/qgis/svg/ User DB Path:
  /usr/bin/qgis/share/qgis/resources/qgis.db

  Provider List Could not find OGR provider! File exists; Path is
  correct ('/home/victorzhiyulee/Desktop/dun.shp', 'dun', 'ogr') Layer
  failed to load!

  Process finished with exit code 0

__author__ = 'victorzhiyulee'
# Importing QGis API
# Importing OGR & OSR
import os
import sys
import PyQt4.QtCore
import PyQt4.QtGui
import qgis.core
import qgis.gui
from qgis.core import *
from qgis.gui import *
from osgeo import ogr, osr
from PyQt4.QtCore import *

# Supply path to the QGis resources on your PC
# noinspection PyTypeChecker
QgsApplication.setPrefixPath("/usr/bin/qgis", True)
# Load providers
QgsApplication.initQgis()
# Show setting of parameters
print QgsApplication.showSettings()

# Load vector layer
data_source = "/home/victorzhiyulee/Desktop/dun.shp"
layer_name = "dun"
provider_name = "ogr"
fileInfo = QFileInfo(data_source)

print('Provider List')
print(QgsProviderRegistry.instance().providerList())

r = QgsProviderRegistry.instance()
if not 'ogr' in r.providerList():
    print 'Could not find OGR provider!'
else:
    print 'Providers found ok!'
# Add layer to the registry
layer = QgsVectorLayer(data_source, fileInfo.fileName(), provider_name)
QgsMapLayerRegistry.instance().addMapLayer(layer)

if fileInfo.exists():
    print("File exists; Path is correct")
    print(data_source, layer_name, provider_name)
    layer = QgsVectorLayer(data_source, fileInfo.fileName(), provider_name)
    if not layer.isValid():
        print("Layer failed to load!")
    else:
        print("Yes, layer loads successfully")
    features = layer.getFeatures()
else:
    print("Check if your path is correct")
    QgsApplication.exitQgis()
iteration = layer.getFeatures()
for features in iteration:
    # Fetch attributes
    attris = features.attributes()
    print(attris)
QgsApplication.exitQgis()


texts: Hide Symbol if the label name not display / rendered in exported file

texts:
I am using QGIS 3.14.
I downloaded world cities name kml, which has markers symbol on every label.
I want to display the symbol only on every name which rendered / displayed on exported file.
Currently, all symbol are displayed even though the cities name are not displayed, so it will only show empty marker without label on it.
How do I filter or hide, so the symbol without names will be hidden in exported file?
As you can see all the symbol which its name not displayed are also displayed on the exported file.


texts: Turn 3D multipolygon with missing faces into cube

texts:
I have data of buildings that are missing walls, like this

MULTIPOLYGON(((5 5 0, 5 0 0, 0 0 0, 0 5 0, 5 5 0),(5 5 5, 5 0 5, 0 0 5, 0 5 5, 5 5 5)))

This would be the floor and roof, stretching from 0 to 5 at a height of 0 and 5 respectively. What I need is to fill the walls in, resulting in the entire object being returned either as TIN or POLYHEDRALSURFACE Z.

I have tried ST_ConvexHull(ST_Collect(geom)) like this

Select ST_asewkt(ST_ConvexHull(ST_Collect(ST_GeomFromEWKT('MULTIPOLYGON(((5 5 0, 5 0 0, 0 0 0, 0 5 0, 5 5 0),
                                                                         (5 5 5, 5 0 5, 0 0 5, 0 5 5, 5 5 5)))'))))

returns only the floor as polygon

POLYGON((0 0 0,0 5 0,5 5 0,5 0 0,0 0 0))

and if I were to use 3DUnion like this

Select ST_asewkt(ST_3DUnion(ST_GeomFromEWKT('POLYGON((5 5 0, 5 0 0, 0 0 0, 0 5 0, 5 5 0))'),
                            ST_GeomFromEWKT('POLYGON((5 5 5, 5 0 5, 0 0 5, 0 5 5, 5 5 5))')
))

it only returns TINs of the floor and the roof

GEOMETRYCOLLECTION(TIN(((5 5 0,5 0 0,0 0 0,5 5 0)),((0 5 0,5 5 0,0 0 0,0 5 0))),TIN(((5 5 5,5 0 5,0 0 5,5 5 5)),((0 5 5,5 5 5,0 0 5,0 5 5))))

Does anyone know how to do this?


texts: Creating a custom OpenStreetMap tile server

texts:
I'd like to create a map tile server based on the entire OSM planet database for a project I'm working on.  I've worked with a couple of guides to try getting something up and running, but haven't found too many comprehensive resources that walk through this process from start to finish.

My requirements so far:

Linux-based (I'd prefer 10.04, but am open to anything)
Tile rendering using Mapnik
The ability to customize the cartography and symbolization to my needs (Quantumnik or Cascadenik seem like good tools for this)

I'd really like to have good custom basemap cartography for some internal apps we're building and I figured this would be a good place to start.  Anyone know of any good guides / resources to building this out?  I'd prefer step-by-step explanatory stuff, as opposed to just the documentation for each of the components...


texts: How to use French LiDAR CSV data in QGIS

texts:
Following on from this Assign elevation data to GPX from LIDAR data
Please could you help with another problem.
I'm trying to use the LiDAR data from France ( http://ids.equipex-geosud.fr/web/guest/france ) in the same manner. But it says it comes with a GeoTIFF DEM file, so it should be easy, but the download doesn't come with one (it does but it's very small file size and QGIS didn't like).
Anyway, it looks like all CSV files. How do I turn them into a DEM? Especially both DTM and DSM?


texts: Multiple value in Query Builder in a Python Script

texts:
I want to do a multiple value selection in this query with an parameter input.
The parameter Regiao_Estado gets a value like the example: [AC, DF] put in the line from the script:

lyr.definitionQuery = 'UF IN' '('Regiao_Estado')'

the whole code is here. The parameter Tabela1 only exists to get the value for the field Ano.
The script will do a query to select values and after that get the sum of the field. For one value it`s ok, But I want for multiple values and maybe multiple fields selection.

Layer1 = arcpy.GetParameterAsText(0)
Tabela1 = arcpy.GetParameterAsText(1)
Regiao_Estado = arcpy.GetParameterAsText(2)
Ano = arcpy.GetParameterAsText(3)
mxd = arcpy.mapping.MapDocument("CURRENT")

for lyr in arcpy.mapping.ListLayers(mxd):
    lyr.name = Layer1
    lyr.definitionQuery = 'UF IN' Regiao_Estado
    arcpy.Statistics_analysis(Layer1, "F:\AgroBD.gdb\Output1", [[Ano, "SUM"]], "")
arcpy.RefreshActiveView()
arcpy.RefreshTOC()

I get the error:

SyntaxError: invalid syntax (M6a.py, line 10)
Failed to execute (Consulta).

I tried many ways to solve this problem but don`t get any progress.


texts: Converting from Ground to Grid coordinate systems in QGIS

texts:
I use QGIS when it's needed, but primarily stick in the AutoCAD world where I have a bit more experience.
I have an orthophoto in Ground coordinates for a project site, and all the survey data that was collected from this stereo orthophoto is ALSO in Ground coordinates.
I wanted to create a CRS in QGIS that is essentially "New Mexico State Plane Central FT, Project" so that I could project other, larger orthophotos and maps that are naturally in "NM State Plane Central FT" into my local coordinate system and have them be accurate locally (i.e., be exportable as GeoTIFF's and importable into my CAD drawings)

Right now, the NMSP data is off by ~250', and I know that my scale factor is 0.999633840 to go from ground to grid, but I can't figure out the appropriate adjustment to my CRS (or an appropriate means to scale a GRID RGB color orthophoto) to get it to line up with my Ground orthophoto.


texts:
Recently, I started making a model on QGIS and now I want to transform it into a plug-in, but it seems that I have to change the code to do it, otherwise it won't work.

What changes are required to perform this conversion?


texts: Satellite image display with the help of GeoServer and OpenLayers

texts:
There is a satellite image. It's size is 10 GB and I need to display this image using GeoServer and OpenLayers. When user select the Satellite image in the layer switcher, I need to display image within 10 seconds.

I tried geopdf but the image quality loss isn't acceptable to customer.

I want to achieve 10 seconds response time using 32 GB satellite image.

Please advice me how to achieve this?


texts: Different color interpreter for Planet Images

texts:
I'm trying to merge a set of Planet Imagery with gdal, but I'm getting Warning 6: gdalbuildvrt does not support heterogeneous band characteristics. Skipping ./out2.tif.

After checking gdalinfo I realized that Blue and Red color are swapped between band 1 and 3 in some images. How is this possible? How can it be fixed?

Actually the solution can be to rewrite all the images using gdal_translate -colorinterp blue,green,red,undefined. But I'm still interested in knowing how this is possible.


texts: Choosing Clip or Insersect tool in QGIS?

texts:
I am doing a analysis but I have information which is not available in all disagregation levels. I mean, some data is at a national level, other is at subnational level and other is in district level. I need all the information in district level.

Can I break down data (national and subnational) into district level (in the atributes table)?

If it is posible with QGIS, should I use clip or intersect?

Weighting with a variable or the size of polygon.

If I do it with clip, atributes table is weighted?

I mean transform data (e.g. In Mexico there are 5,600 diabetes deaths at a national level, I require those data at a subnational level weighted by population, or the size of the polygon).


texts: Creating a point on edge of closest polygon

texts: My project is looking at a grid in the ocean. I want to capture each ocean point's distance to the nearest coastline as a weight, and then use buffer to capture the sum of nightlight data within a 1 degree radius buffer around where the line intersects with the coastline. I know how to work with rasters and buffers but would like to know if there is a way to find the point at which the line crosses the polygon. In the picture below, I have data for the red point and the polygon. I want to find the shortest distance to the polygon, identify where that crosses, and then create a buffer (i.e. the green circle) around that point. The data contained in that buffer would ideally then be reflected in the point.


texts: I have two point layers: warehouses and shops. They have relation like many to many (one warehouse delivers to several shops and one shop connected with several warehouses). I have ID of that connection. I need to show by filter all warehouses, that delivers to one shop or all shops, that connected with warehouse. Now i have a composition of all objects in two layers and using grouping filter in arcgis portal webappbuilder. Is any method to relate two layers for using in group filter without composition? I need to show feature geometry with label, not only table. I'm trying to build relation class, it works fine in desktop, but after publishing on portal there is no related table in map viewer.


texts: Comparing 2 attribute tables to discern unmatched attributes

texts:
In QGIS, I am trying to compare two different attribute tables. One is simply the updated version of the other. Each feature has a unique identifier.

How can I return the number of unmatched records? Is it possible to compare two separate attribute tables in QGIS?


texts:
I'm fairly new to QGIS, but I'm getting the hang of it. I have two layers - one DEM and one raster made from shape file (Buildings), but need one with higher Elevation where the buildings are for terraflow. But I can't seem to get the raster calculator to work with if Statement. I found a Workaround by assigning a fixed elevation value to the building raster, that is higher than any Point on DEM.

But on DEMs with high Elevation difference this means that some Buildings are over 50m tall!!!  would like to make an if statement like this:
 A  B
A = DEM

B = Buildings

if(B>A,A+10,A)


texts: Tabulate intersection of ArcGIS for Desktop returns empty output?

texts:
I am using a US federal lands polygon layer (http://coastalmap.marine.usgs.gov/GISdata/basemaps/boundaries/fedlands/fedlanp020.htm)

I downloaded the layer, reprojected as Albers Equal Area, and then created a 25 km buffer around the polygons. I have been trying to use Tabulate Intersection (ArcGIS 10.2 toolbox), with the buffer as the zones and the original federal lands polygon layer as the classes in order to calculate the percentage area within each buffer that is federal. However, output is empty.

I have tried the following:

Check geometry
Repair geometry
Place both in the same geodatabase
Export as individual shapefiles
subset the buffer feature class (10000 and 1000 rows)
Select all and then run tabulate intersection
Opened new .mxd with only the two input files included

Tabulate Intersection will work when a few records are highlighted, but won't work on any subsets > 1000 records. Overlapping polygons don't seem to be an issue. I know R, so am willing to give that a shot, but have found that it bogs down with this many polygons.

I also attempted to use the intersect tool and received the error messages "ERROR 999999 Error executing function . .Table was not found . . . Invalid Topology". Other stack exchange messages suggest this relates to memory.

Is there a way using Arcpy/Python to iterate over the buffer feature class and perform the calculation on the individual polygons and concatenate the results?


texts: PyQGIS Attribute table show selected features

texts: Is there a way to open attribute table at set the setting to "show selected features" through pyqgis?


texts: Software needed to recreate VBA forms using python for Arc 10.0

texts:
I'm new to programming and need to convert VBA scripted tools for the 10.1 upgrade.

I need to recreate the forms and use Python scripting.

How do I do this?

I need the form to function in 10 so that the upgrade will go smooth.

A list of specific software that I will need would be so helpful.


texts: How to open QGIS attribute table in form view using Python?

texts:
I've got a plugin where I load some layers. After the last layer loaded, I show the attribute table.

However, I would like to select the form view and not the table view (see image).

How can I do this in Python?


texts: How to change symbol size in QGIS map composer legend?

texts:
I need to set an identical size to several symbols inside a QGIS map composer legend. The 'Symbol' menu in the Legend properties allows to set a height and width but this has no effect on the size of the symbols in the legend...

Here is an example of my legend automatically taking symbols sizes from layers:

EDIT : For now, I solved the problem by cloning my layers. Inside each cloned layer, I keep just one record to which I set the symbol size I need in the composer. Then, I build my legend in the composer that uses theses cloned layers. This is just a stopgap solution and I really would like to know if a more direct way does exist with QGIS.


texts: Export selected converted geotagged photo points to a new feature class

texts:
I used the GeoTagged Photos to Points tool to convert a folder of geotagged photos to points, importing the photos into the file gdb as attachments. I now have a feature class of about 7 000 points.

I have selected 100 of the points and now need to put this into a new feature class in a different gdb. However, the points lose their attachments when exported to a different gdb, and copying the original attachment table to the new gdb results in all 7 000 being copied.

Trying to associated the match table with the new feature class using Add Attachments fails because the geotagged conversion tool automatically stores the image paths as relative not absolute. Is there a way for me to extract the selected photos along with their attached images to a new feature class?


texts: How to get shapefiles or kml files for city map in India?

texts:
I am doing a project.

In that I'm using QGIS for visualization purpose.

Using the OpenLayers plug-in i am getting map of a city.

What i want to know is whether I can convert that into shapefiles or KML files?

Or is there any process of getting shapefiles or kml files of a city's roads?


texts: Create legend with OpenLayers 3 / GeoServer

texts:
I am using OpenLayers 3 and GeoServer to create a map. I would like to create a legend for the map that I have created. I've been looking for a good example I could base my map legend off of, but was struggling to do so. I used the following GetLegendGraphic to get the legend I wanted:

http://localhost:8080/geoserver/wms?REQUEST=GetLegendGraphic&VERSION=1.0.0&FORMAT=image/png&WIDTH=20&HEIGHT=20&LAYER=cite:censustractsum&STYLE=basicstyle9&

However, I have no idea how to display this using OpenLayers 3. Is there any way for me to create a legend without having to use anything else. If OL3 can't do this, could I use GeoExt? All the examples of GeoExt I saw were with OpenLayers 2, and I wasn't sure if GeoExt would work with OpenLayers 3.


texts: Inserting space before capital letter in attribute using QGIS field calculator?

texts:
I have a field, and it has many words in it for every row and I want to insert a space before every capital letter except of the first one.

Example:
I have something like this:

  GeorgeJohnMike

or

  George smithJohnMike

And I want to modified it some thing like this:

  George John Mike

or

  George smith John Mike

How can I do it using QGIS or in ArcGis?


texts: I would like to import (in order to create a point layer) in QGIS a csv table with UTMREF/MGRS-coordinates. The coordinates in the column look like this 32UMV6134481745. Usually QGIS requires an X and Y field. Is there a possibility in the import, or do I have to separate or even transfer the coordinates in the csv before import?


texts: How to render a single raster image from defined extent at a specified resolution using raster tiles from openmaptiles as a source?

texts:
I have a problem stitching tiles together to get a specified resolution from map bounds. Im getting the bounds from a leaflet container using openmaptyles raster as a source.

Tile format : (http://host/styles/dark/{z}/{x}/{y}@4x.png), tiles can be sizes of 256,512,1024

Example Bounds : ( Zoom 13, xMin 4148, yMin 2817, xMax 4150, yMax 2819 )

Now, I need to export an image at a resolution of 8268 x 11811px specifically for print.

Ive looked into using mapnik render_to_file to accomplish this but  is it possible to use this source to do it?

The way im doing it now is with a python script which gets the correct tiles from the extent but I cant get the exact resolution im looking for.

I am no GIS expert.


texts: Extract all pixels of given RGB value and render as MultiPolygons

texts:
Given a geotiff image with a colortable in it, what is the easiest way to extract each color as its own "layer or shapefile or multipolygon object"?

Specifically, I am looking for a programmatic approach to solve this problem. (not a GUI app, QGis/ArcGis approach)

Here's what gdalinfo shows:

$ gdalinfo GLOBCOVER_L4_200901_200912_V2.3.color.tif
Driver: GTiff/GeoTIFF
Files: GLOBCOVER_L4_200901_200912_V2.3.color.tif
Size is 129600, 55800
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.001388888888897,90.001388888888883)
Pixel Size = (0.002777777777778,-0.002777777777778)
Metadata:
  Generate by=gdal_mean
  Copyright=Copyright ©UCL Geomatics, BELGIUM 1999-2010
  Authors=Sophie Bontemps <sophie.bontemps@uclouvain.be>,Eric Van Bogaert <eric.vanbogaert@uclouvain.be>,Pierre Defourny <pierre.defourny@uclouvain.be>
  process begin time=2010-12-23T09:49:37
  process finish time=2010-12-23T09:57:38
  process files=CL5_GLOBCOVER-L5_CLASSIF_2009_V2.3.20101220.tif
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0013889,  90.0013889) (180d 0' 5.00"W, 90d 0' 5.00"N)
Lower Left  (-180.0013889, -64.9986111) (180d 0' 5.00"W, 64d59'55.00"S)
Upper Right ( 179.9986111,  90.0013889) (179d59'55.00"E, 90d 0' 5.00"N)
Lower Right ( 179.9986111, -64.9986111) (179d59'55.00"E, 64d59'55.00"S)
Center      (  -0.0013889,  12.5013889) (  0d 0' 5.00"W, 12d30' 5.00"N)
Band 1 Block=129600x1 Type=Byte, ColorInterp=Palette
  NoData Value=0
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
    1: 0,0,0,255
    2: 0,0,0,255
    3: 0,0,0,255
    4: 0,0,0,255
    5: 0,0,0,255
    6: 0,0,0,255
    7: 0,0,0,255
    8: 0,0,0,255
    9: 0,0,0,255
   10: 255,255,100,255
   11: 170,240,240,255
   12: 170,240,240,255
   13: 170,240,240,255
   14: 255,255,100,255
   15: 255,255,100,255
   16: 255,255,100,255
   17: 0,0,0,255
   18: 0,0,0,255
   19: 0,0,0,255
   20: 220,240,100,255
   21: 210,240,100,255
   22: 210,240,100,255
   23: 0,0,0,255
   24: 0,0,0,255
   25: 0,0,0,255
   26: 0,0,0,255
   27: 0,0,0,255
   28: 0,0,0,255
   29: 0,0,0,255
   30: 205,205,102,255
   31: 240,200,100,255
   32: 205,205,102,255
   33: 0,0,0,255
   34: 0,0,0,255
   35: 0,0,0,255
   36: 0,0,0,255
   37: 0,0,0,255
   38: 0,0,0,255
   39: 0,0,0,255
   40: 0,100,0,255
   41: 0,100,0,255
   42: 0,100,0,255
   43: 0,0,0,255
   44: 0,0,0,255
   45: 0,0,0,255
   46: 0,0,0,255
   47: 0,0,0,255
   48: 0,0,0,255
   49: 0,0,0,255
   50: 0,160,0,255
   51: 0,0,0,255
   52: 0,0,0,255
   53: 0,0,0,255
   54: 0,0,0,255
   55: 0,0,0,255
   56: 0,0,0,255
   57: 0,0,0,255
   58: 0,0,0,255
   59: 0,0,0,255
   60: 170,200,0,255
   61: 0,0,0,255
   62: 0,0,0,255
   63: 0,0,0,255
   64: 0,0,0,255
   65: 0,0,0,255
   66: 0,0,0,255
   67: 0,0,0,255
   68: 0,0,0,255
   69: 0,0,0,255
   70: 0,60,0,255
   71: 0,0,0,255
   72: 0,0,0,255
   73: 0,0,0,255
   74: 0,0,0,255
   75: 0,0,0,255
   76: 0,0,0,255
   77: 0,0,0,255
   78: 0,0,0,255
   79: 0,0,0,255
   80: 70,75,0,255
   81: 0,0,0,255
   82: 0,0,0,255
   83: 0,0,0,255
   84: 0,0,0,255
   85: 0,0,0,255
   86: 0,0,0,255
   87: 0,0,0,255
   88: 0,0,0,255
   89: 0,0,0,255
   90: 40,100,0,255
   91: 40,100,0,255
   92: 40,100,0,255
   93: 0,0,0,255


texts: Explain attribute selection for arcpy.SpatialJoin_analysis

texts:
I have exported a ModelBuilder model to python script, but it is failing.  I suspect it is the attribute selection that is the problem.  This model does a series of Spatial Join tasks one after the other (5 in total).

I have created a simple script that does a single Spatial Join to use as a test.  This works successfully, but I do not understand what the various inputs mean.

Can anyone explain what the options are that I don't understand.

# Import arcpy module
import arcpy

# Script arguments
Property_SpatialJoin2 = arcpy.GetParameterAsText(0)
if Property_SpatialJoin2 == '#' or not Property_SpatialJoin2:
  Property_SpatialJoin2 = "filepath\\Property_SpatialJoin2" # provide a default value if unspecified

# Local variables:
Property_shp = "filepath\\Property.shp"
MM_EIL_shp = "filepath\\MM_EIL.shp"

# Process: Spatial Join
arcpy.SpatialJoin_analysis(Property_shp, MM_EIL_shp, Property_SpatialJoin2, "JOIN_ONE_TO_ONE", "KEEP_ALL", "HAS_Code \"HAS_Code\" true true false 30 Text 0 0 ,First,#,filepath\\Property.shp,HAS_Code,-1,-1;Contact \"Contact\" true true false 50 Text 0 0 ,First,#,filepath\\Property.shp,Contact,-1,-1;EIL_vers \"EIL_vers\" true true false 10 Long 0 10 ,First,#,filepath\\MM_EIL.shp,EIL_vers,-1,-1", "INTERSECT", "", "")

The section that I do not understand is what the options are for the attributes to include:

HAS_Code \"HAS_Code\" true true false 30 Text 0 0 ,First,#,filepath\\Property.shp,HAS_Code,-1,-1

The sections I (think I) know are:
 First refers to the join match, and can also be Min/Sum/Mean etc.
 HAS_Code \"HAS_Code\" gives the name of the attribute as it is called in the existing layer followed by what will be used in the new layer
 filepath\\Property.shp,HAS_Code gives the path to the file that has the attribute

I am not sure what the true true false 30 Text 0 0 refers to
or what the # and -1,-1 do.


texts: Distance between two lat/lng-points (geography) on a plane

texts:
I use PostGIS with PostgreSQL to calculate the distance between two places, specified by a latitude and longitude coordinate. The outcome of 18.5 km is what I expect:

SELECT ST_DISTANCE(epicenter, house_location)/1000 "distance_in_km"
FROM (
  SELECT ST_SetSRID(ST_Point(53.345, 6.672),4326)::geography(point,4326)         epicenter
  ,      ST_SetSRID(ST_Point(53.2081669, 6.5759146),4326)::geography(point,4326) house_location
) x;

What I want to know is the distance between those two points, without considering the curve of the earth.

My use case: I have a location of an earthquake (epicentre) and a location of a house. I also know the depth of the earthquake, and with Pythagoras I calculate the distance between the house and the hypocenter.

These calculations are Euclidean, but the distance calculation is not. I've heard the difference is minimal, but I want to know for myself what the difference is, how I can calculate it?

I assume both places have the same distance to the earth's center (but if you know how to take that fact into the equation, I would be interested as well).


texts: Split or hide specific data points from CSV file in QGIS

texts:
I have about 100 locations that are open for a couple of weeks, some for a day or two, some for multiple days, some every day. They are currently all plotted as data points on a map. So right now, I have a map with every location that will be open on every day over the weeks.
What I need to do is separate these points so that I can create specific maps for each day. For example, I need a map of all locations that will be open on Day 1, all locations that will open on Day 2, and so on.

I realize that I can go into my CSV file that I used to plot these points and manually separate them into multiple files based on each individual day (or build an ungodly Excel macro). However, since I am trying to increase my QGIS skills, I am wondering if there is a more efficient way that I can solve this problem within QGIS itself.

(I understand this may be possible by using the Python plugin, but I am very new to Python and cannot really do anything with it, as I have only started to pick it up in the last few days. So if you were to recommend using Python, please do so only if you have a script)


texts: Convert high resolution raster (tiff) file into single polygon?

texts:
My objective is to clip a single shapefile based on multiple tiff files. Consider that we have a map of a country and we want to clip the states that have residential areas. The residential areas are represented by multiple tiff files.

My current method employs converting tiff into shapefile using the Raster to Polygon tool and using dissolve tool to get a single polygon (iteratively done with multiple tiffs and merging them together to get the final polygon). However this route takes a long processing time and results in freezing/ sudden exit from ArcGIS 10.4.1 application.

Current procedure

Raster-Polygon : Raster is single banded with values ranging from 0-250 with 30m x 30 resolution. The resulting polygon is very large with several polygons.
Dissolve (each tiff)
Merge
Clip the single shapefile using the merged polygon feature

I am assuming there should be an efficient way to capture the boundary/ extent into a polygon and then merge these polygons, instead of converting and dissolving the polygons from each raster file.


texts: Adding multiple Excel files and display XY data using ModelBuilder

texts: I would like to use ModelBuilder to add multiple Excel files and display the XY data as a shapefile and eventually convert the point data to raster. I am using ArcGIS Desktop 10.4. I have tried the "Make XY Event Layer" only one layer displays.


texts: ogr2ogr - conversion shapefile to kml - how to reduce the size of the output file?

texts:
Is possible to reduce the size of the output.kml file?

ogr2ogr -f "KML" output.kml input.shx


texts: Loading data using QGIS to PostgreSQL database

texts:
I am trying to import a shapefile using desktop QGIS to postgreSQL database using the import shape file to postgresql database tool but it raises the following error:

  Problem inserting features from file:
  E:/MOA/sample point identification/EthioWoredaandZones_csa2007/EthioWoredas_csa2007.shp

The database gave an error while executing this SQL:

INSERT INTO "public"."EthioWoredas_csa2007"("REGIONNAME","ZONENAME","WOREDANO_","WOREDANAME","the_geom") VALUES ('Addis Ababa','Region 14','140104','B... (rest of SQL trimmed)

The error was:

  ERROR:  Geometry type (Polygon) does not match column type (MultiPolygon)

How can I import this shapefile into postgeSQL and manipulate it?


texts: Remove first character from a string if it satisfies an if statement

texts:
I have a bunch of phone numbers and some contain a 1 in front of them e.g. "18885556607" and I would like to simply remove the 1. I can't just remove all first characters because there are numbers like "8665558411"

Here is the code I've tried:

def splitme(s):
    if (s[:1] == "1"):
        return s[2:]
    else:
        return(s)

expression:

splitme( !Phone! )

The field is a string field, please let me know what I need to change.


texts: Building raster (population density) dependend diagram in QGIS

texts:
I'm using QGIS and trying to solve the following problem.
I got 2 layers:

population density raster
point layer with the information about number of people could be hosted/served etc in this points.

I want to show the area where such number of people (that could be
hosted/served) lives around this points (for instance, using circle diagrams).
I thought about creating multibuffer, then counting the number of people living
in each of them and manually trying to find the most appropriate size of it. However, I believe there is more rational option.
Could you please propose any solutions?


texts: Error when installing the "Minkowski fractal dimension calculator" Plugin in QGIS 3

texts: I have tried to install the Minkowski fractal dimension calculator via the QGIS Plugin Repository, but got this error shown in the screenshot.


texts: Cartographic technique for detailing a 'Focus Area'

texts:
I'm working on updating a planning document.  Whoever was in charge of cartography last time around used a great technique to identify and showcase focus areas (where to focus money and effort for maximum return on investment).  I'm trying to use this same type of technique by interpreting and recreating the work-flow.

As you can see, it certainly tells more of the story than a standard flat-filled polygon especially when considering the landscape from a restoration perspective.

The focus area seems to be a terrain / colored hillshade -- could I recreate this with the national elevation dataset (NED) from USGS?  Suggestions for a solid workflow to get me going would be great.


texts: LeafletJS: Load Tile Layer based on bounding box

texts:
I have a LeafletJS map application with lots of layers. I have a few layers that only have imagery for certain regions across the world. Because everything goes on a world base map. All tile caches are generated as standard tile caches covering the world in web Mercator (Google tile cache system) and are served using a ArcGIS Server.

Here is what I would like to do: Load a tile layer when certain bounding box comes into active view and then unload it when box is out.

Is there a simple way to do this in LeafletJS without going through manual geometry check ?


texts: Error from aggregate function from the R package sf

texts:
I am currently reading the new (and free!) book by Edzer Pebesma and Roger Bivand (you can find it here).

In chapter four they describe how one can aggregate all the measurement stations within one state. I follow there code but get an error. How so?

library(sf)
library(spacetime)
library(rgeos)
library(stars)
# dataset from spacetime package
data(air)
d = st_dimensions(station = st_as_sfc(stations), time = dates)
aq = st_as_stars(list(PM10 = air), dimensions = d)
a = aggregate(aq, st_as_sf(DE_NUTS1), mean, na.rm = TRUE)

The last line returns the error:

Error in aggregate.data.frame(as.data.frame(x), ...) :
Arguments must have same length


texts: Geoserver "too many clients already" error

texts:
I'm creating an application to display some data with geoserver and postgres database.
My application is working fine when I display 5-6 layers of data one by one, after that, the next layers I try to display are not displayed correctly, only a few tiles are on the screen. When I look up the Geoserver database I have error message

  FATAL: sorry too many clients already connected


texts: Finding and using mean and standard deviation value in field calculator

texts: I am creating a field for a z-score of a particular attribute. This requires calculating the mean and the standard deviation for the attribute across all features. I've been using the statistics window to get the mean and standard deviation, copying them into a field, and then doing the calculations in a third field, but I feel like that's not very efficient. I can't find a way to calculate the mean using the field calculator and my scripting skills are pretty weak.


texts:
Working with OpenLayers 2.11, I'm desperately trying  to change my feature style as the user selects it (fairly simple I agree!).

Indeed, It seems easy but I can't see where It fails, perhaps someone experienced with OpenLayers could enlighten me.

Below is my code :

   function onFeatureSelect(event)
   {

    var feature = event.feature;

    var styleWaypointSelect = new OpenLayers.StyleMap(
    {
        externalGraphic: 'img/map_marker.png',
        graphicWidth: 16, graphicHeight: 26, graphicYOffset: -24,
        title: 'selectedFeature'
    })

    var clonedFeature = waypointLayer.getFeatureById(feature.id);
    clonedFeature.style = styleWaypointSelect;

    waypointLayer.drawFeature(clonedFeature,styleWaypointSelect);
    waypointLayer.redraw(); /* refresh() I tried several things to refresh the layer but it doesn't seem to change anything at all*/

EDIT :

I finally figured out that OpenLayers StyleMap class allows different properties, "default", "delete", "temporary",  "select".
(you can see this when you debug the style object in firebug)

I added the "select" style of my vector layer like this :

        var styleWaypoint = new OpenLayers.StyleMap({
              "default": new OpenLayers.Style({
                externalGraphic: 'img/mobile_loc.png',
                graphicWidth: 16, graphicHeight: 26, graphicYOffset: -24,
                title: 'marker'
            }),
               "select": new OpenLayers.Style({
                externalGraphic: 'img/mobile_loc_select.png',
                graphicWidth: 16, graphicHeight: 26, graphicYOffset: -24,
                title: 'marker_select'
            })
        });


texts: Save style/renderer in the QGIS project file

texts:
I have a vector layer in QGIS and using pyqgis I have created two styles/renderers:

style_1 = QgsFeatureRendererV2(...)
style_2 = QgsFeatureRendererV2(...)

vector_layer.setRendererV2(style_1)

One of the styles is assigned to the vector layer and the other one is for later use.

Saving the QGIS project will store style_1 with the layer but the second style will be lost if the project is closed. I would like to store the second style in the QGIS project file, not in a seperate style file.

Using QgsProject.writeEntry(...) I can store numbers, texts and lists in the QGIS project file.

Question:
Is there a corresponding method for storing custom styles/renderers in the QGIS project file?


texts:
I have two large datasets, RP1 and RP2. I would like to subtract one from the other, however the datasets are extremely large and so exit out with a MemoryError when I try a direct read:

RP1 = rasterio.open(big_file_path)
RP1 = RP1.read(1, masked=True)

RP1 = RP1_dataset.read(1, masked=True)

  File "rasterio_io.pyx", line 348, in

  rasterio._io.DatasetReaderBase.read MemoryError

Therefore I have tried to resort to windowed reading and writing as explained by the rasterio docs. I wrote the following script:

import numpy
import rasterio

image = numpy.ones((150, 250), dtype=rasterio.ubyte) * 127
RP1_filename = r'RP1.tif'
RP2_filename = r'RP2.tif'

RP2 = rasterio.open(RP2_filename)
profile = RP2.profile
with rasterio.open(RP1_filename) as RP1:
    for block_index, window in RP1.block_windows(1):
        RP1_block = RP1.read(window=window,masked = True)
        RP2_block = RP2.read(window=window, masked = True)
        result_block = RP2_block - RP1_block
        print(result_block.shape)

RP2.close()

Since RP1 and RP2 have the same window schema I took the liberty of indexing into RP2 by the current RP1 block index. This seems to work fine and the result_block variable contains a masked array as expected.

The problem is I don't know how to collect all these arrays properly such that I can write a new output raster with the result values.

Does anyone know how to do that?


texts: Calling gdalogr:rasterize algorithm in python console in QGIS?

texts:
I'm trying to rasterize a simple line-vector layer. I can do it by using QGIS GUI, but I want to do the same with python console.
Unfortunately it always gives back an "Error: Wrong number of parameters" or "Wrong parameter:___".

Call like processing.alglist('gdalogr:rasterize') gives the shown list of parameters, where for many of them there is no further explanation what they are and how to set them. I'm using QGIS 2.16.1.

Has anyone made a successful call for this algorithm?

ALGORITHM: Rasterize (vector to raster)

INPUT <ParameterVector>
FIELD <parameters from INPUT>
DIMENSIONS <ParameterSelection>
WIDTH <ParameterNumber>
HEIGHT <ParameterNumber>
RAST_EXT <ParameterExtent>
TFW <ParameterBoolean>
RTYPE <ParameterSelection>
NO_DATA <ParameterString>
COMPRESS <ParameterSelection>
JPEGCOMPRESSION <ParameterNumber>
ZLEVEL <ParameterNumber>
PREDICTOR <ParameterNumber>
TILED <ParameterBoolean>
BIGTIFF <ParameterSelection>
EXTRA <ParameterString>
OUTPUT <OutputRaster>


texts: Watersheds and Catchment delineation in Earth Engine only with included datasets

texts:
Is there a possibility to derive Watershed boundaries or catchment areas by using one or more of the integrated data sets in google earth engine?

The available data sets are:

WWF HydroSHEDS Hydrologically Conditioned DEM, 3 Arc-Seconds
WWF HydroSHEDS Drainage Direction, 3 Arc-Seconds
WWF HydroSHEDS Void-Filled DEM, 3 Arc-Seconds
WWF HydroSHEDS Flow Accumulation, 15 Arc-Seconds

I know i can export them as raster and work with QGIS or ARCMap to derive the boundaries but anybody knows how to derive it with python or java script inside the google earth engine without exporting the data?


texts:
I have a raster file (TIF file) and want to use its extent as a mask. I need this mask to be able to extract other shape files /raster files by this mask.

How can I convert my raster file to a simple mask (shp.file) regarding only the outline of my raster file?


texts: Geodatabase to PostGIS character encoding

texts:
Similar to another question I posted previously, I'm trying to import datasets from Esri geodatabases into a PostGIS database directly, and inside the data there are greek letters and symbols that are not converting correctly.  I tried using FWTools on Windows with a command like this:

ogr2ogr -f "PostgreSQL" PG:"host=server.domain.com user=gisuser dbname=gis"
dataset.mdb LayerName -a_srs EPSG:4326

The first time I tried, I got this error from ogr2ogr:

ERROR 1: INSERT command for new feature failed.
ERROR: invalid byte sequence for encoding "UTF8": 0xb4
HINT: This error can also happen if the byte sequence does not match the
 encoding expected by the server, which is controlled by "client_encoding".

I found out that setting a particular environment variable at the Windows command line can fix this issue:

SET PGCLIENTENCODING=LATIN1

That fixed importing problem -- all the records and geometry are in PostGIS, but as a result, the non-english characters are not encoded correctly.  When I use UTF-8 or UTF8 for client encoding (not sure which is correct), I get the same invalid byte sequence error as above.  Has anyone else had this problem or found a way around it?  Is this an issue with the .mdb file itself, or is there something on the server I can change to globally correct this?


texts: Adding point shapefile to map shuts down program

texts: I am trying to add a point file that I have used before to a map.  Every time I add it the program (ArcGIS) crashes.  I can add other files, including the same data for a different year, but this specific file crashes the program, so I assume there is something wrong with the file. It would be very difficult to recreate the file, so I need to diagnose what might be wrong and then repair the file...if possible.


texts: "DBF file can not be opened" error using shp2pgsql

texts:
I am trying to import a shapefile using shp2pgsql. I wrote this command:

shp2pgsql -I -s 2154 -W "LATIN1" W:\IGN\BDPARCELLAIRE-RGF93-LB93-MAJ2013-V4\BDPV_1-2_SHP_LAMB93_D01\COMMUNE.shp CADASTRE.COMMUNE | psql -U postgres -d omexom

But, when I execute this, I have an error :

Unable to open W:\IGN\BDPARCELLAIRE-RGF93-LB93-MAJ2013-V4\BDPV_1-2_SHP_LAMB93_D01\COMMUNE.shp or W:\IGN\BDPARCELLAIRE-RGF93-LB93-MAJ2013-V4\BDPV_1-2_SHP_LAMB93_D01\COMMUNE.SHP.
W:\IGN\BDPARCELLAIRE-RGF93-LB93-MAJ2013-V4\BDPV_1-2_SHP_LAMB93_D01\COMMUNE.shp: dbf file (.dbf) can not be opened.

I really don't understand, because I tried to import the same layer using the graphical interface PostGIS and it works.

Where am I wrong in my code?

The dbf file is present.


texts:
I've been looking for a manner to change a vector extent in QGIS, but it's unsuccessfull!

I want change extent because in the past, I had problems to convert vector to raster. I try present the problem bellow:

I'm InVEST user. InVEST is  a set of models developed by Natural Capital Project to mapping, quantifying and valuing Ecosystem Services. A while ago, I had a problem with the Carbon model. Actual LULC and Future LULC are input datas beside others. When I analysed this datas I perceived different extents. The model uses the extent when converting vector to raster during processing, therefore the Actual LULC pixels and Future LULC pixels did not fit properly. Another model's process is map algebra. When model executing map algebra, the results present a pattern due to LULC datas misalignment. I corrected this problem using the same extent for two LULC data.In ArcGIS the way to acomplish this task is very simple.

I was ArcGIS user. I've been using QGIS recently. ArcGIS and QGIS work with extent in different ways? Are there a manner to change vector extent in QGIS?


texts: ArcObjects: How to get the size of the formatted text part only of a TextElement whose text contains formatting tags?

texts:
I'm working in ArcGIS 10.0 and Python 2.6.5.  I have a string, say 'mMx', that I want to enter into a TextElement (in the data view).  The lowercase letters are to be substituted with a graphical symbol through use of the usual html-style formatting tags.  We'll assume that the regular font is Arial (specified in the TextSymbol) and that the sybolic font is Wingdings.  So, I build a formatted string like so:

<FNT name="Wingdings" size="48">m</FNT>M<FNT name="Wingdings" size="48">x</FNT>

This string is what I set as text for my TextElement.  Now, what I want to get is the bounds of the formatted text only, as it would appear on the screen with full symbology.  The text appears fine -- with symbology -- on the screen.  The problem is that the bounds of the text element include the length of the whole input string, formatting tags and all.  In map units, the text is maybe 50' long, but the bounds are 1000'.  I have tried QueryBounds and QueryOutline (both of which say they take symbology into account) and also ITextSymbol.GetTextSize.  They all return the same crazy dimensions.  ArcGIS is processing the formatting tags just fine as far as the text itself is concerned, but the bounds include the tags.  I don't get it.  This is a big concern, because I absolutely need to know the onscreen-text-only size of the TextElement for use in subsequent calculations.

Does anyone have any insights on this problem?  Thank you all for the help ... and I apologize for asking such a specific ArcObjects question here, but I get exactly NO response from ESRI forums anymore.

Solution found:

Well, I just found the problem, and it has to do with the way the format string I show above was made in Python.  I have been doing a lot of HTML code block output in Python lately, and I have been used to using triple quoted strings.  That is what I did when building the format string.  So the following is the actual line of code that builds the format string.  (Lesson learned.  I should always include a code snippet.)

fmtchar = """<FNT name="Wingdings" size="48">""" + eachChar + """</FNT>"""

It should be:

fmtchar = '<FNT name="Wingdings" size="48">' + eachChar + '</FNT>'

I find it interesing that ArcGIS would process the formatting tags and also count the tag characters as part of the text.  Hmmm.  BTW, let me know if I should put this solution in an answer to my own question.


texts: Joining table to feature class using ArcGIS with Delimiter?

texts:
I have a table of attributes that correspond to different cities in my state in an excel file that I am trying to join to a polygon feature class, based on the city name.

After executing the join, when I go to symbolize based on these new attributes, this is what I see under the symbology tab:

I cannot list the unique categorical values; it looks like a delimiter issue where several columns are somehow being combined, but there are no extraneous commas, punctuation marks, etc. in the table.  The table originated in Excel; I saved as a .csv then used the 'Table to Table' tool to import it into ArcGIS as a table in a file geodatabase.

Where in the join process did I go wrong?


texts: How to use Facebook (GeoTIFF) maps on Carto?

texts:
My intuition is that the maps that Facebook generated for most countries are most useful on interactive maps like Carto's, as the killer feature is the resolution for which you'd need to zoom in.

However, I am not familiar with GeoTIFF, let alone importing them to Carto.

What would be an example workflow if someone would like to see, say, children under five of Hungary overlayed on a basemap?


texts: Send Geometry to query using ArcGIS REST API

texts:
I am writing a windows console app in VS 2010 to perform an overlay using the ArcGIS Rest API.

My first REST query returns information about the geometry for a specific feature.  I then want to send that geometry to another feature class, to return data about the shared area.

I am using the DeserializeObject method to convert the query results to a Dictionary.  How do I take the geometry from that result set and post it to a new query?

Update:

What if I user the JavaScriptSerializer.Serialize method on the geometry string?  Is there a way to pass the resulting JSON string to the REST service call?  I already tried just assigning it to the "geometry" value in the querystring, but that failed.


texts: QGIS/QField - Joining updated and new features to original shapefile

texts:
I am currently working in QGIS 2.18 with Qfield to collect address information.

I currently have a list of addresses (.csv) that were given as part of the project which contains a variety of information, including point XY co-ordinate for each address (although these appear to be wrong and need geocoding.)I have made this into a shapefile for use in QField offline editing.

Samplers will use Qfield and walk around each town, moving the existing address point features to within buildings polygons at the correct addresses and updating the attributes of each point with the required information. They will also be adding new point features for those addresses which are currently not in the existing address list.

Each evening, I download the QField shapefiles from each of the samplers containing the updated existing points and the new points that they have taken in the day.

If it was just a case of updating the old points with the new attribute information, I understand that I can join across a common field and repopulate each field with the new information provided. I have tried this, however with multiple files to add (5 samplers/folders a day) and manually repopulating the old columns with the updated fields is taking an enormous amount of time especially as I need the new XY coordinate for each feature.
It also means that my original shapefile is getting confusing as I have so many joined columns.

As they will also be adding new point features, I cannot find a way to add these new points without updating the .csv itself and re-adding it to QGIS as they do not have a common column with the original .csv to be included in the join above.

Can anyone suggest the ideal/most efficient way of using the shapefiles from each sampler to update the attributes for each existing point feature as well as including the newly created ones?

I chose QField as QGIS is the software i'm most comfortable working in. If there is a significantly easier alternative to Qfield for gathering and integrating the new data then yes I am more than happy to try others but they would have to be free.


texts: Problems in producing animation GIF files in R

texts:
I am trying to produce a GIF file for my maps in R using gganimate but errors show up.

devtools::install_github("dgrtwo/gganimate")

library(ggplot2)
data(sp6, package = 'aqp')
p <- ggplot(sp6, aes(x = pH, frame = hz)) +geom_density()
library(gganimate)

Error in fetch(key) : lazy-load database 'C:/R/R-3.4.2/library/gganimate/help/gganimate.rdb' is corrupt

(I think this can be fixed by restart R session? Anyway, the real problem below ...)

gganimate(p)

Executing:
"C:\ImageMagick-7.0.7-Q16\convert.exe -loop 0 -delay 100 Rplot1.png Rplot2.png Rplot3.png
    Rplot4.png "file148c3014687f.gif""
'C:\ImageMagick-7.0.7-Q16\convert.exe' is not recognized as an internal or external command,
operable program or batch file.
an error occurred in the conversion... see Notes in ?im.convert
Error in file(file, "rb") : cannot open the connection
In addition: Warning messages:
1: running command 'C:\WINDOWS\system32\cmd.exe /c "C:\ImageMagick-7.0.7-Q16\convert.exe" --version' had status 1
2: running command 'C:\WINDOWS\system32\cmd.exe /c "C:\ImageMagick-7.0.7-Q16\convert.exe -loop 0  -delay 100 Rplot1.png Rplot2.png Rplot3.png Rplot4.png "file148c3014687f.gif""' had status 1
3: In cmd.fun(convert) :
  '"C:\ImageMagick-7.0.7-Q16\convert.exe -loop 0  -delay 100 Rplot1.png Rplot2.png Rplot3.png Rplot4.png "file148c3014687f.gif""' execution failed with error code 1
4: running command '"C:\ImageMagick-7.0.7-Q16\convert.exe -loop 0  -delay 100 Rplot1.png Rplot2.png Rplot3.png Rplot4.png "file148c3014687f.gif""' had status 127
5: In file(file, "rb") :
  cannot open file 'C:\Users\gjiang\AppData\Local\Temp\RtmpA95rFH/gganimate\file148c3014687f.gif': No such file or directory

Is there a way to fix this?

P.S

After installing the software ImageMagick 7.0.7, I used p <- animate(ras.stack) in the raster package and then used saveGIF(p) in the animation package, but still getting

'C:\ImageMagick-7.0.7-Q16\convert.exe' is not recognized as an internal or external command, operable program or batch file.


texts: Deleting the old layer set automatically - after sorting layers of one specific group within the table of contents alphabetically

texts:
I'd like to reorder the layers within ONE specific group of the TOC alphabetically by layer name, but without reordering the toplevel layers or the layers of other groups.
Gérman Carillo, thanks again, helped me to solve this (see link below)
I'm now struggling with removing old layers automatically (I just have some basic knowledge in R, but none in Python):
In the sample of Gérman Carillo the old layer structure is automatically removed, and only the layers in the intended order remain (see animation on earlier question of 27.5.2021.). This is exactly the result I'd like to have.But if I run the script in my QGis (3.16), always both remain, one set of all layers in the origin order and another set in the new/intended order (see image below).Anybody an idea what's going wrong here? Or how to modify the script to delete the former layer set?
This is the current code:
from collections import OrderedDict

root = QgsProject.instance().layerTreeRoot()
group = root.findGroup("group1") # name of the to-be-sorted-group
reverse_order = False

LayerNamesEnumDict=lambda listCh:{listCh[q[0]].name()+str(q[0]):q[1]
                                   for q in enumerate(listCh)}

mLNED = LayerNamesEnumDict(group.children())
mLNEDkeys = OrderedDict(sorted(LayerNamesEnumDict(group.children()).items(), reverse=reverse_order)).keys()

mLNEDsorted = [mLNED[k].clone() for k in mLNEDkeys]
group.insertChildNodes(0,mLNEDsorted)
for n in mLNED.values():
    group.removeChildNode(n)

This is the image of the result after running the script


texts: Google Earth Engine: Finding least cloudy pixel of Landsat 8 images in ROI

texts:
I am having trouble cloud-screening Landsat 8 imagery in Google Earth Engine for a country with high cloud cover throughout the year.

I have tried various techniques but unfortunately there is a substantial amount of cloud in all imagery and the techniques given in the Earth Engine Help are not quite cutting it.

I've tried simply requesting the median of each band as given here but clouds remain. The Google Earth Engine SimpleCloudScore script works much better but leaves some gaps. The code that has worked the best for me is the code given by Nicholas Clinton in this answer. Here it is in full:

    var ic = ee.ImageCollection("LANDSAT/LC08/C01/T1_RT_TOA");

    // A polygon representing the roi.
    var geometry = ee.Geometry.Polygon(
            [[[-121.85897778617999, 37.70881514186375],
              [-121.83975284337708, 37.76202899390253],
              [-121.94137041100134, 37.759857750255144]]]);

    var c = ic.filterBounds(geometry);

    var withCloudiness = c.map(function(image) {
      var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
      var cloudiness = cloud.reduceRegion({
        reducer: 'mean',
        geometry: geometry,
        scale: 30,
      });
      return image.set(cloudiness);
    });

var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));
print(filteredCollection);

This issue is that this only brings up entire images for a small ROI or single point of interest.

What would be ideal is if this code could be edited so that it ran the same calculation per pixel, selecting the most cloud-free image for each individual pixel within the time range given. I've tried combining it with the SimpleCloudScore script so that a function calls back to the extended simpleCloudScore algorithm in Richard's code but no luck so far.

Of course, if anybody has any completely different approaches to suggest it would be great to hear them also.

I was hoping this question could produce useful answers to all new users building up their scripts for cloud-screening in Earth Engine.


texts: ArcGIS 9.3. for Silverlight Missing color swatch inside legend

texts:
I have an ArcGIS 9.3. MapService and a Silverlight Application that consumes it. The application is deployed on two diffrent Systems. On system A everything is fine. On system B the map service is displayed correctly, as well as the legend for it. But what is missing are the color swatches inside the legend. I can see that the color swatch images are created inside the arcgisoutput folder. And I can access them via the web browser. When I debug the application of system B with fiddler, I see that there are no request trying to get the swatch images, as it is the case for system A.

What could be the problem?


texts: GeoPandas gdf.geometry.buffer for spatial join

texts:
In that part of the code I'm creating geodataframes (only point features) out of two CSV files. I want to merge all points from the two dataframes within a range of 150m. Therefore I create a buffer and then do a spatial join. However I need the original df later. Therefore I make a copy of the db before buffering. The join is working, but afterwards the original df (hikr_gdf) is containing polygons, not points. I don't get it, since I am not working on hikr_gdf, but on buffer_gdf? Probably I'm overseeing something simple, but what could be wrong?

files = ["../spatial_data/hikr_waypoints.csv", "../spatial_data/osm_waypoints.csv"]

gdfs = []
for file in files:
    df = pd.read_csv(file, delimiter=",")

    gdf = GeoDataFrame(
        df.drop(['latitude', 'longitude'], axis=1),
        crs={'init': 'epsg:4326'},  # wgs84
        geometry=[Point(xy) for xy in zip(df.longitude, df.latitude)])
    gdfs.append(gdf)

    gdf.to_file(f"{file[:-4]}_processed")

hikr_gdf = gdfs[0]
osm_gdf = gdfs[1]
hikr_gdf = hikr_gdf.to_crs({'init': 'epsg:3857'})
osm_gdf = osm_gdf.to_crs({'init': 'epsg:3857'})
hikr_gdf.to_file("../spatial_data/hikr_gdv.shp")

buffer_gdf = hikr_gdf
print(f"print1:\n{hikr_gdf.geometry}")

buffer_gdf.geometry = buffer_gdf.geometry.buffer(150)
print(f"print2:\n{hikr_gdf.geometry}")

hikr_gdf.to_file("../spatial_data/test.shp")

waypoints_gdf = gpd.sjoin(osm_gdf, buffer_gdf, how='left', op='within', lsuffix="osm", rsuffix='hikr')

The problem is between the two print-statements. The output is:

print1:
0      POINT (1238304.930 5936603.562)
1      POINT (1240420.095 5937074.189)
2      POINT (1241967.945 5936348.037)
3      POINT (1241604.869 5936252.491)
4      POINT (1258256.795 5973883.395)
                    ...
478    POINT (1426811.788 6067915.272)
479    POINT (1430399.552 6066195.439)
480    POINT (1436366.417 6008753.012)
481    POINT (1437629.379 6017117.804)
482    POINT (1436299.730 6016980.603)
Name: geometry, Length: 483, dtype: geometry
print2:
0      POLYGON ((1238454.930 5936603.562, 1238454.208...
1      POLYGON ((1240570.095 5937074.189, 1240569.373...
2      POLYGON ((1242117.945 5936348.037, 1242117.223...
3      POLYGON ((1241754.869 5936252.491, 1241754.147...
4      POLYGON ((1258406.795 5973883.395, 1258406.073...
                             ...
478    POLYGON ((1426961.788 6067915.272, 1426961.066...
479    POLYGON ((1430549.552 6066195.439, 1430548.829...
480    POLYGON ((1436516.417 6008753.012, 1436515.695...
481    POLYGON ((1437779.379 6017117.804, 1437778.657...
482    POLYGON ((1436449.730 6016980.603, 1436449.008...
Name: geometry, Length: 483, dtype: geometry

I don't get why hikr_gdf is changing, when I am buffering buffer_gdf.


texts: Subtract two DEM to find volumetric change for a glacier without bedrock topography

texts:
I am a masters student, and trying to calculate the volume change for a glacier. I have 4 different DEMs, none of them are of the bedrock topography below the glacier. I have read that you can simply find the volume change by subtracting DEMs (older-newer) and then multiply by cell size to get the volume change. But the numbers I am getting does not give much sense to me at this moment, it basically tells me that the glacier is gaining volume, rather than losing it.

Anyone have any experience doing this type of analysis?


texts: Convert large ASCII file to an image file

texts:
I have a file of 5gb in a .asc form. The first lines are

ncols         43200
 nrows         18000
 xllcorner     -180
 yllcorner     -60
 cellsize      0.0083333333333333
 NODATA_value  -9999

It is a global layer.
Trying to convert it to .tif is not possible since maybe it is to big.
I converted it using QGIS and raster --> translate to a .nc file, but it seems that the result is not as expected. I get values  for the north hemisphere above 20Lat only but not for the south(all value have the same.

Any suggestions on how can I convert it and in what form?

In the properties in Qgis the file has the following information:

CRS
EPSG:4326 - WGS 84
Geographic
Extent
-180.0000000000000000,-60.0000000000000000 : 180.0000000000000000,90.0000000000000000
Unit
degrees
Width
43200
Height
18000
Data type
Float64 - Sixty four bit floating point
GDAL Driver Description
AAIGrid
GDAL Driver Metadata
Arc/Info ASCII Grid
Dataset Description
C:/lbii.asc
Compression

Band 1

More information

Dimensions
X: 43200 Y: 18000 Bands: 1
Origin
-180,90
Pixel Size
0.008333333333333333218,-0.008333333333333333218

When converting it to a .img i get the following image(almost same as .nc)


texts:
I have Linux Mint 14 KDE which is based on Ubuntu 12.10 Quantal.
In the QGIS download page, there's no "Quantal" repositiries. I'm trying to use the "Precise" repositories, but I'm not allowed to install the latest QGIS version (1.8) and I received many messages of broken packages.

Any suggestion?

Thanks!

Lucas


texts: Disconnect users who are creating locks on specific feature classes

texts:
I have a nightly Python script which updates some domains in a geodatabase. The domain owner is "AMIS". In addition to the AMIS schema in the geodatabase, there are other schemas. In order to update the domains, I must disconnect any users which are connected to data which may be using these domains. The data which may be using these domains are only in the AMIS schema.
I would like to disconnect any users which may be connected to the AMIS data. It's okay if they're also connected to other datasets.
In the example below, you can see that the lock owner connecting to AMIS data are the users "AMIS" and "MANAGEMT".

I would like to disconnect these two users, but not any other users, such as "ASSET", "VISITOR", "COMMERCIAL".
Currently, I'm using arcpy.ListUsers to get a list of all users that may be connected to the database. However, I may inadvertently disconnect a user that isn't locking any "AMIS" data.
I'm also aware of the arcpy.TestSChemaLock, but that only returns a boolean. I would like to know who is locking the data so I can disconnect them.
How can I determine which users are locking data which begins with "AMIS", that is, data in the AMIS Schema?
I'm using ArcGIS 10.8. Geodatabase is Oracle backend.


texts: Using Field Calculator to Convert Data Type In Field

texts:
I am trying to convert a field that is listed as "String" because it was imported into QGIS as a CSV (rural poverty population percentages by county in NC). I need to convert the percentages field from "string" to "Real" with at least one tenths decimal place. I have been using the field calculator to create a new field with a data type conversion (toreal( "Field_Name")) with the precision at 2 (for 2 decimal places I am assuming). The expression is valid and I am able to get a result as a new field. The problem is that I keep getting resulting values from the adjacent column to the one that I am trying to convert. I get real numbers with 2 decimal places but from a column next to the one I am trying to convert! (This adjacent field happens to be a "string" data type (RCU codes) so I am definitely not mistaking the data in it's converted state.

Has anyone had a problem with the field calculator and it selecting a different field that the one you have identified in the conversion expression?


texts: Create a slide to control opacity for layers

texts: I am starting with my programming skills. I m self learning so I am stuck in  this point I am trying to create a web map where charge few different layers (wms) but I would like to give functionality to them providing control with the  opacity so you can overlay two layers at same time, but I have no idea how to do it. I tried few ideas I read from here but I am not sure what I am doing wrong.


texts:
Is there anyway to style base maps using Leaflet?  I know there are quite a few options for loading different pre-made styles (such as Stamen, Mapbox, etc...), but I am looking for a way to create my own style without going back to Google Map's API.  Or is there a way to style one of the particular base maps?

Basically, I'd like to keep using Leaflet for this project as it is quite nice and simple, but without that ability I'll have to move on.  Thanks.

Edit:  Is there a way, through Leaflet, to use Google Maps custom styling?  I like the coding aspects of Leaflet as compared to GMaps, but still would like the flexibility in terms of design that custom styling of the tiles gives me.  As this is a very small project, I don't want to set up my own Tile Map Server if possible.


texts: Projection between EPSG:4326 and EPSG:4269 using pyproj does nothing

texts:
According to this quite thorough answer :
Is re-projection needed from SRID 4326 (WGS 84) to SRID 4269 (NAD 83)?
the crs EPSG:4326 and EPSG:4269 are slightly different, due to a difference in one of the semi-axis of the reference ellipsoids (small difference) and also to a shift in the origin of the NAD83 Datum used in 4269.
However, when using pyproj (python 3.8.8 and pyproj 3.0.1 on windows) to project a point from EPSG:4326 to EPSG:4269, I obtain the exact same lat/lon :
from pyproj import CRS, Transformer
from shapely.geometry import Point
from shapely.ops import transform

crs_4326 = CRS("EPSG:4326")
crs_4269 = CRS("EPSG:4269")

transformer_4326_to_4269 = Transformer.from_crs(crs_from=crs_4326, crs_to=crs_4269).transform
p = Point(-118,33)  # some position in California
p1 = transform(transformer_4326_to_4269, p)

print(p)  # returns POINT (-118 33)
print(p1) # also returns POINT (-118 33) !

I am quite confused as I was expecting a (rather small) difference.


texts:
I am trying to convert my old ArcGIS 9.3 (VBA) code to ArcGIS 10 using Visual Basic 2008 Express Edition.  So far I have sucessfully created an pop-up window that displays the attributes of a selected feature in text boxes on the form.

I added a reference to the editor and made some adjustments to my code. My form still crashes when I try to press the button. Am I missing something or am I just referencing the attribute table in correctly?

THANKS!!!

Private Sub cmdUpdate_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles cmdUpdate.Click
    pMxDoc = My.ArcMap.Document
    Dim pMap As IMap
    pMap = pMxDoc.FocusMap
    Dim pFLayer As IFeatureSelection
    Dim Count As Integer

    'Get a reference to the editor.
    Dim editorUID As ESRI.ArcGIS.esriSystem.UID
    editorUID = New ESRI.ArcGIS.esriSystem.UID
    editorUID.Value = "esriEditor.Editor"
    Dim editor As IEditor
    editor = My.ArcMap.Application.FindExtensionByCLSID(editorUID)

    'Get FeatureLayer
    pFLayer = pMap.Layer(Count)
    Dim pFCursor As IFeatureCursor = Nothing
    pFLayer.SelectionSet.Search(Nothing, False, pFCursor)
    Dim pF As IFeature
    pF = pFCursor.NextFeature

    'Get Dataset from FeatureLayer
    Dim pEditDataset As IDataset
    pEditDataset = pFLayer

    'If not editing, start editing dataset's Workspace
    If editor.EditState <> esriEditState.esriStateEditing Then
        editor.StartEditing(pEditDataset.Workspace)
    End If

    'Initialize first (only) feature
    Dim pFeat As IFeature
    pFeat = pFLayer.NextFeature

    'Update attributes with user input
    pFeat.Value(pFeat.Fields.FindField("YEAR")) = txtPrjYear.Text
    pFeat.Value(pFeat.Fields.FindField("PROJECTYPE")) = txtPrjType.Text
    pFeat.Value(pFeat.Fields.FindField("LOCATION")) = txtLocation.Text
    pFeat.Value(pFeat.Fields.FindField("AGENCY")) = txtAgency.Text
    pFeat.Value(pFeat.Fields.FindField("DESCRIPTION")) = txtDescrip.Text
    pFeat.Value(pFeat.Fields.FindField("COST")) = txtCost.Text
    pFeat.Value(pFeat.Fields.FindField("DATEEDITED")) = txtDate.Text
    pFeat.Value(pFeat.Fields.FindField("USEREDITED")) = txtUser.Text

    'Stop editing and save edits
    editor.StopEditing(True)

    'Refresh map
    pMxDoc.ActiveView.Refresh()

    'Close form
    Me.Close()
End Sub


texts: Get the length of all lines in Python QGIS?

texts:
How can I get the length of all lines with python in QGIS.

For example if I have:

How can I say, I have 23 m (example) of lines?

I found that if I click in one part of the line:

But really I wouldn't know how can I get the Length(Cartesian) or Length(Ellipsoidal) in Python.

Because I do the next code but I don't find any fields with the name "length" but in the above image I can see that inside of "Derived" There is a field.

print(features[0])
print(features[0].geometry())
lista = features[0].fields().toList()

And for example when I right click in the layer, and I open "open attribute table" I see left image, but on the other hand I show that the same line has a field "length" how you can see in the right image:


texts:
Is there any way of doing the "keep only matching records" on QGIS when working with joint tables?

(This is an option on ArcGIS, but not on QGIS as far as I know)


texts:
I want to use open street map as cascading WMS to my Geoserver. Which capabilities URL should i use?

Can any body tell me the capabilities URL for OSM ?

Regards


texts: Creating polygons with "Create shapefile layer" always produce multipolygons instead of polygons when saving in QGIS

texts:
Not sure if is a bug or I'm missing something.
Every time i'm creating a new shapefile with "Create shapefile layer" i got a multipolygon instead of a polygon (at least is what I'm expecting).
Also got same issues,  when using an import of PostGIS table (with "Polygon" type geometry in db) and trying to add/save the polygons, I got an error with the message of unable to commit because different type inside the database compared to the default type which is created.
I ve tried different versions of QGIS (3.10,3.12,3.16) but seems to have all buggy behaviour.
I'm on Windows.

Creating a simple polygon as a Polygon type with rectangle tool

Before saving: Getting the type of the polygon displays correctly a Polygon type

After saving: The type is changing in Multipolygon without any warning.


texts: Geoserver CQL Cross Layer Intersection Failing

texts:
I am having trouble forming a properly formatted CQL request using geoServer cross layer filer functions and was hoping others could assist. My request is as follows:

http://localhost:8080/geoserver/wfs?request=GetFeature&typeName=AU:TownPoints&CQL_FILTER=INTERSECTS(the_geom,collectGeometries(queryCollection('AU:Districts','the_geom','CODE=A1')))&version=1.1.0

The intent is to select the towns that are within a selected District shape.

With the above request a null data set is returned:

<wfs:FeatureCollection xmlns:ogc="http://www.opengis.net/ogc" xmlns:AU="hidden.com" xmlns:tiger="http://www.census.gov" xmlns:cite="http://www.opengeospatial.net/cite" xmlns:nurc="http://www.nurc.nato.int" xmlns:sde="http://geoserver.sf.net" xmlns:wfs="http://www.opengis.net/wfs" xmlns:topp="http://www.openplans.org/topp" xmlns:it.geosolutions="http://www.geo-solutions.it" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:sf="http://www.openplans.org/spearfish" xmlns:ows="http://www.opengis.net/ows" xmlns:gml="http://www.opengis.net/gml" xmlns:xlink="http://www.w3.org/1999/xlink" **numberOfFeatures="0"** timeStamp="2013-03-20T22:36:01.593Z" xsi:schemaLocation="hidden.com http://localhost:8080/geoserver/wfs?service=WFS&version=1.1.0&request=DescribeFeatureType&typeName=DHS%3ATownPoints http://www.opengis.net/wfs http://localhost:8080/geoserver/schemas/wfs/1.1.0/wfs.xsd">
<gml:featureMembers/>
</wfs:FeatureCollection>

I think the formatting of my query string 'CODE=A1' is wrong, so geoServer returns a null set for the queryCollection, meaning the intersection is being performed against a null set. (Note: CODE='A1' works fine if I query Districts directly.)

I have also tried 'CODE='A1'' and some other variations but this nesting always causes an exception similar to the following:

Could not parse CQL filter list. Encountered "A1" at line 1, column 92.
 Was expecting one of:
 ")" ...
 "," ...
 "/" ...
 "*" ...
 "+" ...
 "-" ...
 Parsing : INTERSECTS(the_geom,collectGeometries(queryCollection('AU:Districts','the_geom','CODE='A1''))).

I have looked at lots of examples but cannot find one with a filter string in the queryCollection function.


texts: Where can I find documentation on QGIS C++ Plugin Development?

texts:
Is external C++ plug-in support available for QGIS?

I'm having a really hard time getting started here. I followed the doc here: http://download.osgeo.org/qgis/doc/manual/qgis-1.1.0_coding-compilation_guide_en.pdf

But that's fairly outdated. I was able to compile and link, but the plugin isn't being brought into QGIS - so I'm wondering if that tutorial is out of date.

When I go to the wiki page at QGIS's website - they reference a VERY outdated page that is identified as obsolete: http://hub.qgis.org/wiki/quantum-gis/Writing_C++_Plugins

I can't seem to find a current tutorial that describes C++ plugin support for QGIS 1.8 (or any version close to it). I realize Python is more pervasive, but at this point we have specific needs that are supported via C++.


texts: How to activate a FeatureLayer operation from a button using dojo/on function?

texts:
this is the setup of the FeatureLayer:

 on(MallaFeatureLayer.getFeatureLayer(), "before-apply-edits", lang.hitch(this, this.asignarUnidadLecturaMalla));

how could I link it to a button?


texts: String value as a number in modelbuilder

texts:
I have got a drop-down list of values in my model. There are word combinations that do not exist in my table.

One of the tools in my model should select string data of a field by number of characters (Char_length). My string variable is called like "Roads type". So I need to have something like

CHAR_LENGTH (" VEHICLES ") = %Roads type%

So can I make a match of words in drop-down list and specific number only within a model without making of additional layers?


texts: Converting LinestringM to Polygon PostGIS

texts:
I'm working with PostGIS and plotting polygons on a Google Earth plugin. I've found a few examples via Google on how to convert a LinestringM to a Polygon for rendering, but none of them, including the PostGIS docs, define very well the params for a LinestringM. Can someone identify what's happening here? Specifically, I need to know what each of the four params of LinestringM are supposed to represent.

Snippet:

....
parameters.put("wkb_geometry", "ST_MakePolygon(ST_GeomFromText('LINESTRINGM(75.15 29.53,77 29,77.6 29.5, 75.15 29.53)'))");
....


texts:
I have two rasters, r1 and r2. I have numerous NA in both of them. I would like to set NA in r2, where r1 NA occurs. Meaning for every NA that occurs in r1, I would like to set NA to r2.

I want to, set the values of r2  to NA that are also NA in r1.

r1 #example raster
     1  2  3  4  5  6  7  8  9   10
1    3  4  NA 6 NA 8 NA NA NA NA NA
2    NA NA NA NA 9 NA 8 NA NA NA NA
3    11 NA NA 18 12 NA NA NA NA  NA

r2
  1  2  3  4  5  6  7  8  9 10
1 3  4  4 6 NA 8 NA NA NA NA NA
2 NA NA NA 19 20 8 NA NA NA NA
3 11 41 NA 18 12 NA NA 23 NA  NA

I am trying this code, but I am pretty sure this code is different of what I want to achieve. Help would be greatly appreciated

 values(r2)[!is.na(values(r1))] <- NA


texts: Is there a way to do an Intersect in SQL, that not only returns a geom column, but also all the attribute information?

texts:
I am trying to do an intersect in SQL using the programme pgAdmin 4. I have managed to do an intersect which only returns the resulting geom column.

I would like to have all the Attribute Information from all the other columns, too though. Just like when you do an intersect in QGIS or ArcGIS.

Here is what I have used that got me only the geom column (I am using colours to make it less complicated and easier to read):

CREATE TABLE public.green AS
SELECT ST_INTERSECTION (b.geom, y.geom_2)
FROM public.blue b, public.yellow y
WHERE ST_INTERSECTS(b.geom, y.geom_2) = 't';

I know the public isn't necessary here, but I have more than one Schema that I use. So I thought it would be better to add the Schema here so all parts are accounted for.


texts: Add more than one layer of GeoJSON data to Leaflet?

texts:
This is my code for loading one layer of .geojson data to Leaflet. The data is basically in the form of markers and it's demonstrating different location in Saudi Arabia. It's saved as "saudi_location.geojson". And my code succeeds in showing this layer on the map.

Say I have another .geojson file called "saudi_airports.geojson" and I want to load it to the map too... and I want the user to be able to toggle between the two layers, maybe check/uncheck a box to show/hide the airports layer or the location layer.

How can this be achieved?

<!DOCTYPE html>

    <title>Leaflet</title>

    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.0.3/dist/leaflet.css" />


texts: Negative values in area calculation as Geopackage

texts:
I'm getting strange results with area calculations via the field calculator, when I'm working with a geopackage file.  The field calculator is returning a negative values for some of the polygons.  I can get it to return the (apparently) correct values by saving the layer as a shapefile, and using the same calculation.

Is there a way to correct this behavior without exporting the layer to a different format, or do I need to avoid geopackages when manipulating data?

This is similar to the problem in some other posts:
Why do area calculations return negative values?
Area is calculating wrong using $area in field calculator (QGIS 2.8.1-Wien)
Wrong $area from field calculator (QGIS)

The file I'm having trouble with is a selection of parcels from a larger parcel data set with parcels within a limited extent merged by owner:
https://www.dropbox.com/s/nwgg8stxrglcdn8/MergedParcels.gpkg?dl=0

There are multi-part features in the file, but I still get the negative values after using multi-part to single part vector geometry tool.

I'm using QGIS 2.18.2 on a Windows computer. My project is in EPSG: 102748 (Washington State Plane North- US Feet), as is the geopackage.

I'm getting the same results with On-the-fly projection both on and off.
The formula I'm using in the field calculator is:

round($area/43560,2)

My project properties are set to calculate areas and distances in map units (feet/square feet).


texts: Using QGIS field calculator with geometry operations between layers

texts: A curiosity question. I am looking for making a field calculation with QGIS, I see is possible to include the geometry field of the same feature layer, nevertheless I was wondering if it is possible to use the geometry from other layers features in my calculation. For example, intersection or else.


texts: Making icons and text larger on high resolution screen for QGIS?

texts:
I have a new laptop with a very high resolution screen and can't get QGIS to work properly.

Icons and text are too small.

I tried to change the resolution, but had no luck.

I also tried to change settings in QGIS but that changes only text and some of the icons.


texts: Listing all ArcGIS Coordinate Systems using ArcObjects?

texts:
We currently have an ArcMap addin that will convert selected or all point features to a SEGP1 file. We also allow projections to a different Coordinate System and this is where it gets tedious. We want to list out all valid coordinate system from arcgis. Currently, we have an xml file where we store the Coordinate System list.

Is there any way where we can retrieve the list of coordinate systems using ArcObjects?

We are using ArcObjects 10.1 .NET btw.


texts: Calculate raster pixel maximum accross multiple bands

texts:
I am using QGIS 3.10 and try to solve the following task. I have a multiband raster (4 bands in total), where each band represents the spatial extent of a tropical storm and the respective pixel value represents the wind speed.
In other words, the multiband raster contains all tropical storms within a given area for a respective year. In this case there have been 4 tropical storms, represented by the 4 bands. Now I want to create a singleband raster containing the max value of the 4 rasters, so that the new singleband raster image contains the maximum wind speed exposure for cell i in year t. If a cell was hit by two storms, I want to assign the value of the band with the higher wind speed.
The question in short is: How do I tell QGIS to compare the raster values of all bands and then choose the highest one. I tried the raster calculator, but I could not figure out a suitable option and I am not really familiar with the syntax in the raster calcaaltor.


texts: QGIS problem exporting to both SVG and PDF

texts:
I'm having a problem when exporting to SVG and PDF from qGIS. When I open the SVG or PDF in Illustrator, the path of my map is severely cut off with lots of data missing (see below).

When I open the PDF in Acrobat, it looks fine. So this is a problem with how Illustrator is handling the files.

I've tried simplifying the geometry and can't get it down below ~70,000 vertices. I've also tried splitting the object in two but then the exported SVG is completely blank.

Does anybody know how to solve this problem?


texts: QGIS Web Client in Windows does not load map

texts:
I have installed QGIS (1.8 Lisboa) desktop and server on a Windows XP machine. I am now trying to configure the QGIS Webclient. I downloaded the zip file from GitHUb, and made the changes to the templates, etc. When I open the qgiswebclient.html page the map does not load. The toolbars and other parts of the GUI are all there.

From OSgeo4w Shell I checked my mapserver is running: Mapserver v5.6.6

I copied the qgis-web-client folder into the osgeo4w\Apache\htdocs folder.

Does anyone know which parameters I should check?


texts: Export single band from a multiband stacked raster with rasterio

texts:
I have been trying to break down my multi-band raster into individual bands with rasterio and I am getting this error:

ValueError: Source shape (11356, 9488) is inconsistent with given indexes 1

here's the code:

import os
import rasterio

 src = rasterio.open(wv_tif)
    for band in range(1, src.count):
        single_band = src.read(band)

        # get the output name
        out_name = os.path.basename(wv_tif)
        file, ext = os.path.splitext(out_name)
        name = file + "_" + "B" + str(band) + ".tif"
        out_img = os.path.join(folder, name)

        print(out_img + " done")

        # Copy the metadata
        out_meta = src.meta.copy()

        out_meta.update({"count": 1})

        # save the clipped raster to disk
        with rasterio.open(out_img, "w", **out_meta) as dest:
            dest.write(single_band)

I have no idea what is going wrong


texts: Seeking example .PRJ file using WKT syntax for Concat_MT?

texts:
I've been digging into the WKT format definition found here:WKT Format with the aim of defining a PRJ with two chained affine transforms to go from a known projected system to a local grid.

I've been unable to generate any kind of PRJ that works with packages like GDAL and was wondering if anyone has any kind of example file that uses the CONCAT_MT entity alongside a standard projection?

Even getting any working implementation of Concat_MT with GDAL at all would be a start, I've tried a few variations like the one below with no success.

PROJCS[
    "WGS84 like system to local ground coordinates",
    GEOGCS[
        "Unknown datum based upon the WGS 84 ellipsoid",
        DATUM[
            "Not_specified_based_on_WGS_84_ellipsoid",
            SPHEROID[
                "WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]
            ],
            AUTHORITY["EPSG","6030"]
        ],
        PRIMEM[
            "Greenwich",0,
            AUTHORITY["EPSG","8901"]
        ],
        UNIT[
            "degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]
        ],
        AUTHORITY["EPSG","4030"]
    ],
    CONCAT_MT[
      PARAM_MT["Mercator_2SP",
        PARAMETER["semi_major",6370997.0],
        PARAMETER["semi_minor",6370997.0],
        PARAMETER["central_meridian",180.0],
        PARAMETER["false_easting",-500000.0],
        PARAMETER["false_northing",-1000000.0],
        PARAMETER["standard parallel 1",60.0]],
        PARAM_MT["Affine",
        PARAMETER["num_row",3],
        PARAMETER["num_col",3],
        PARAMETER["elt_0_1",1],
        PARAMETER["elt_0_2",2],
        PARAMETER["elt 1 2",3]]]
    UNIT["metre",1,AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]
]

Does anyone have a working example PRJ with multiple affine transformations applied to a projected system?


texts:
I'm reading a very large json file as an InputStream in Android with JsonReader. This json has a geometry object like so..

 "geometry":{
      "type":"MultiPolygon",
      "coordinates":[
         [
            [
               [
                 -7.661981835965626,
                 39.1261981217978
               ],
             [
                 -7.661983712214173,
                  39.12617044304507
             ],

           ]
        ],
(...)

What's the best way of parsing this geometry into a WKT string?

I'm using something like this for the other fields, but for the geometry it's a bit complex and confusing.

if (tag.equals("geometry")) {
    reader.beginObject();
    while (reader.hasNext()) {  (...)


texts: Can vector grid shapefile be clipped with satellite image in QGIS?

texts:
I am working in OpenLayers plugin in QGIS. In order to project the satellite image with respect to vector shapefile I have, Andre Joost suggested me to create a vector grid with CRS system of vector shapefile. As suggested I created grid.

Now my problem is this:

Am I able to clip the vector grid shapefile with the satellite image obtained?


texts: Geoserver security and application security separation

texts:
I am using geoserver that publises WMS, WFS services to clients. I am using geoserver services in an OpenLayers application. My application will requiere username and password authentication.

1 - For example I have a polyline data published WMS. If a user is working on Arizona, she/he can see only Arizona polylines. Can I configure this on geoserver. Or should I do on application?

2- Geoserver admin panel is including Users/Groups/Roles menu. Should I use this users and roles in OpenLayers applicaiton?

How can I specify the security way?


texts: How to port QGIS C++ plugins to other QGIS installations?

texts:
I developed a C++ plugin for QGIS using VC++ express 2008 and a QGIS dev release (1.9.90 alpha), my plugin works fine on this release and now I want to port it to other QGIS installations, so I copied the generated dll to the plugins directory of a QGIS 1.7.4 installation but it didn't get recognized by the extension manager.

Anybody can help with this?

Thx!


texts: QGIS crashes when I try to open the extension tab

texts:
As I try to boot up the manager in order to manage my currently installed plugins, QGIS tries to load some of my custom data/plugins repositories from work. Those need some specific credentials to access, they work on other computers.

As QGIS tries to load the first one, it freezes as it asks for my IDs. What I am looking for is a simple fix without needing to delete and install QGIS back again, maybe by somehow erasing the repositories without opening the extension tab or trying to load them for instance.


texts: Batch Processing and Raster Calculators in QGIS

texts:
As of right now I have two sets of raster files. The first and second set are of differing sizes. I am attempting to use set to clip the other, using the process outlined in this question here Clipping raster layer with *raster* mask layer in QGIS (the one that says you can "create a unitary mask on the fly ")
This process works. However, the QGIS raster calculator is unable to use batch processing, and my understanding of the GDAL calculator is it can only be used on similarly sized raster files. I have attempted to use the SAGA raster calculator, but I am unable to get any kind of equation that works.
What equation would I need to use in the SAGA raster calculator to produce similar results and, if this is not possible, what other QGIS specific options do I have?


texts: Converting 3D kml files to 3D Shapefiles using ArcMap without 3D Analyst

texts:
How do I convert 3D KML files to 3D Shapefiles using ArcMap without 3D Analyst?
I currently do not have access to 3D analyst. I am using ArcGIS Desktop 10.4.1.
We are migrating from Google Earth EC to ArcGIS Portal


texts: How can I split my polygon file into topological boundaries?

texts:
This question was asked 8 years ago, but the algorithm in the accepted answer is gone. (How do I split polygons into line segments?)
So how can I split my polygon file into "segements"? The lines should be broken up like this picture shows:


texts:
I store my geographical information in PostgreSQL database.

I want to use GraphHopper for my routing propose, but GraphHopper can read information in shapefiles or OSM format.
I now convert my information to shapefile format and use GraphHopper.
But I want to read information directly from PostgreSQL and for this propose I have to write a custom reader for this.

How can write this reader?


texts: Aggregate by date-time and compute a centroid as an output?

texts:
(I would prefer using R if possible)

I have a table with coordinates (in latitude, longitude), time information (YYYY-mm-dd HH:MM:SS) and ID information (example provided). I sometimes have multiple coordinates in a one minute window for one ID. I would like to produce a table with continuous time (by="min") with:

NA: if no coordinates are recorded at this particular time this particular ID.
coordinate: if there is only one coordinate in the one minute window.
centroid: if there are multiple coordinates

Example coordinate table:

Coordinates<-data.frame(GMT_Date= c("2017-04-16 15:15:00","2017-04-16 15:15:00","2017-04-16 15:15:00","2017-04-16 16:33:00","2017-04-16 16:33:00","2017-04-16 16:33:00","2017-04-16 16:33:00"),
                        Latitude= c(15.0,15.2,15.4,16.0,16.5,16.3,16.2),
                        Longitude= c(30.2,30.1,30.2,31.5,31.6,31.8,31.9),
                        ID= c("id1","id1","id1","id2","id2","id1","id1"))

Example continuous time table:

startDate<-as.POSIXct("2017-04-16 00:01:00",format="%Y-%m-%d %H:%M")
endDate<-as.POSIXct("2017-05-23 23:59:00", format="%Y-%m-%d %H:%M")
time.table<-data.frame(GMT_date= seq(startDate,endDate,by="min"))
head(time.table)

Expected final output:

output.expected<-data.frame(GMT_date=c("2017-04-16 00:01:00","2017-04-16 15:15:00","2017-04-16 16:33:00"),
id1_centroid= c("NA","Latitude, Longitude", "Latitude, Longitude"),
id2_centroid=c("NA","NA","Latitude, Longitude"))


texts: Add headers to QGIS symbology and add colour ramp

texts:
I have added a catagorical symbology based on a simple filter in QGIS Desktop 3.6.3.
I have three different columns on which the coloring of features should be based. Combinations of column values are possible and not all columns have data.
The symbology is based on the following expression, I am content with the output:
CASE

    WHEN Veg_codom <> '' AND  Veg_onderg <> '' THEN Veg_dom + '+' +         Veg_codom
    WHEN  Veg_codom IS NULL AND  Veg_onderg <> '' THEN Veg_dom
    WHEN Veg_codom <> '' AND  Veg_onderg IS NULL THEN Veg_dom + '+' +       Veg_codom
    WHEN Veg_codom IS NULL AND Veg_onderg IS NULL THEN Veg_dom

END

The layout now looks like this:

I am struggeling for a while now to produce the following, but no succes:

So: I want to add headers to my symbology (that also translate to my legend for print-composer) and I want to be able to control colour ramps for each group. In Arcmap this is fairly easy to do and this is actually the first thing I find harder to find out on QGIS.
Atrribute table is as follows:


texts:
I'm trying to generate the legend of my mapguide map from the mapguide API but I've found any samples...

Have you an idea of how to use GenerateLegendImage function on a layer and how to convert the MgByteReader result to base64 img string ?


texts: Which GIS software or plugins can perform a Cellular Automata Simulation?

texts: I am interested in running cellular Automata to run a simulation on city growth. Which GIS software can run this out of the box, or via plugins?


texts: How to create a grid to do kriging interpolation that correspond to areas in a map having actual observations?

texts:
I need to do kriging interpolation but wanted to have the interpolation in areas where I have actual data; not across the whole map.
This is my example.
Here is the base map (shp file in a .zip file): http://www.filedropper.com/lpr000b16ae
Here is the data I want to interpolate: http://www.filedropper.com/yielddata
library(dplyr)
library(ggplot2)
library(sf)
library(gstat)
library(stars)

## a) Defining the projection of the maps
proj_lambert <- "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

## b) Map Canada selected provinces.
canada.map <- st_read("C:/Users/rotundjo/Desktop/lpr_000b16a_e/lpr_000b16a_e.shp") %>%
  filter(PRNAME == "Manitoba" |
           PRNAME == "Alberta" |
           PRNAME == "Saskatchewan") %>%
  st_transform(proj_lambert )

## c) Data set
data.yield <- read_csv("C:/Users/rotundjo/Desktop/yield.data.csv")
data.yield.coords <- data.yield %>% select(lat, long)
data.yield.data <- data.yield %>% select(av.yield)

#### Creating LAT LONG SpatialPoints
coordinates(data.yield.coords) = c("long", "lat")
proj4string(data.yield.coords) <- CRS("+init=epsg:4326") # 4326 is the EPSG identifier of WGS84.

### Converting the LAT LONG to the lambert CRS shown above
lambert_data.yield.coords = spTransform(data.yield.coords, proj_lambert)
lambert_data.yield.coords = as.data.frame(lambert_data.yield.coords)

#### Adding the Type column back to the data frame, with the new polar coordinates
lambert_data = cbind(lambert_data.yield.coords, data.yield.data)

coordinates(data.yield) = c("long", "lat")
proj4string(data.yield) <- CRS("+init=epsg:4326")
data.yield = spTransform(data.yield, proj_lambert)
data.yield.sf <- st_as_sf(data.yield, coords = c("long", "lat"), crs = proj_lambert)

## d) Variogram fit
v = variogram(av.yield~1, data.yield.sf)
plot(v, plot.numbers = TRUE)

v.m = fit.variogram(v, vgm("Sph", psill = 3, range = 3000000, nugget = 1))
plot(v, v.m, plot.numbers = FALSE)

## e) Building a empty grid
st_bbox(canada.map) %>%
  st_as_stars(dx = 1000) %>%
  st_set_crs(proj_lambert) %>%
  st_crop(canada.map) -> grd

## f) Kriging interpolation
k = krige(av.yield~1, data.yield.sf, grd, v.m)

## g) plotting all together
library(viridis)
ggplot() +
  geom_stars(data = k, aes(fill = var1.pred, x = x, y = y)) +
  geom_sf(data = canada.map, color = "gray", fill="transparent") +
  geom_point(data = lambert_data, aes(x = long, y = lat), size = 0.1, color = "gray") +
  theme_bw()

library(RColorBrewer)

ggplot() +
  geom_stars(data = k,
             aes(fill = cut(var1.pred,breaks=c(0,20,21,22,23,24,Inf),
                            labels = c("Below 20",
                                       ">20-21",
                                       ">21-22",
                                       ">22-23",
                                       ">23-24",
                                       "Above 24")),
                            x = x, y = y)) +
  scale_fill_brewer(palette = "Spectral") +
  geom_sf(data = canada.map, color = "gray", fill="transparent") +
  geom_point(data = lambert_data, aes(x = long, y = lat), size = 0.1, color = "gray") +
  theme_bw()

This is my result:

But I would like to have the interpolation made only in the areas having actual observations. That is, the grd that is defined above should be defined by an area defined by the actual points. The expected result should be something like this:

The white stripped are should be part of the grd for the interpolation.
I was not able to define the grid area for interpolation to be restricted somehow to the areas where there are actual observations.


texts: Difference between X,Y and Latitude,Longitude or Meters

texts:
According to Esri:

"To add a table of x,y coordinates to your map, the table must contain
two fields: one for the x-coordinate and one for the y-coordinate. The
values in the fields may represent any coordinate system and units
such as latitude and longitude or meters."

I am getting confused here.
For example, on GPS what is X,Y and are they different than Latitude,Longitude?


texts: Querybuilder in QGIS-->select features in a shp-file, which touch others

texts: I have a ESRI-shp-file called veg.shp and the field BYKP. Is it possible to make a query within the querybuilder and select all polygons with the beginning attriubte "F" and the second condition where they touch all polygons with the beginning attriubte "O" in the field "BYKP"


texts: Making Open Street Maps Transparent in Leaflet

texts:
Info:

I have a local Open Street Maps Server which is serving tiles for Michigan, Indiana, Ohio, Illinois, and West Virginia. I also have a geoserver set up serving a wms layer containing states. Because I only have several states served in Open Street Maps, OSM does not show state borders, which is why I have a WMS layer.

The Problem:

I wanted to layer OSM over top of my WMS layer and make it transparent so that the state lines will show through from the WMS layer.

Code:

var map;
var feature;

function load_map() {

     map = new L.Map('map', {
            zoomControl: true,
            center: [0, 0],
            zoom: 3
     });

    var osmUrl = 'http://my_server/osm_tiles/{z}/{x}/{y}.png',
    osmAttribution = 'Map data &copy; 2012 OpenStreetMap contributors',
    osm = new L.TileLayer(osmUrl, {
            maxZoom: 18,
            attribution: osmAttribution,
            transparent: true
    });

    var states = L.tileLayer.wms('http://my_server:8080/geoserver/Indigital/wms', {
            format: 'image/png',
            transparent: false,
            layers: 'Indigital:states'
    });

    var overlays = {
            "States:WMS": states,
            "OSM": osm,
    };

    var layerControl = L.control.layers(null, overlays);
    map.addControl(layerControl);