texts: What is the purpose of GeometryCollection Option Associated with the Shape Dropdown?

texts: I have some data sets that have point, line and polygons in the same table and I have been selecting the Type "GeometryCollection" for the shape column.  I have noticed that items such as labels do not behave consistently when I make this selection.  What is the purpose of this option and is there a better way to handle multiple Shape types?


texts: Can't insert a value into the table. PyQGIS 3.0. Postgresql

texts:
Can't insert a value into the table. PyQGIS 3.0

import psycopg2

conn = psycopg2.connect(dbname='postgres', user='postgres', password='1', host='localhost')
cursor = conn.cursor()

cursor.execute(""" INSERT INTO table_name (name, make, model, year)
                    VALUES('Alex', 'Toyota', 'Supra', '2017-10-12');
                """)

cursor.close()
conn.close()

I run the code, but nothing data is added.


texts:
I need to model a fault plane below the earth's surface and find distances from different points on the surface of the earth.

I know the dimensions of the plane and its orientation but I can do that on a paper.

How to set this up in ArcGIS and find the distances I need?


texts: QGIS form control to show multiple selection data

texts:
I need to show labels for coded values recorded in one single field, i.e.:

Crops Field
Row 1: 1 2 3
Row 2: 1 3
Row 3: 1 4

Codes:
1 - Maize
2 - Beans
3 - Sugar Cane
4 - Others

I want to show the labels for this data within a QGIS form, when a user clicks for information, is this possible, I can´t find an appropriate controller.


texts: Extend polygon to polyline feature (keeping attributes)

texts:
I am looking to extend the polygon feature to the to the coastline polyline (blue line) while keeping the attributes of the polygon features (Figure below).

What I am looking to do is find the percent ownership (public vs. private lands) along per coastline mile (i.e. I only care about the right-most parcels of the polygon. And I do have these separated out). A simple buffer does not retain attributes for me to know the ownership of the adjacent parcel to the left of the coastline. Also since the coastline meanders an intersect does not capture all the parcels.

In one attempt, I converted the polygon to a raster, so that I could use the "Expand" tool to extend out to the line. But when I did so, the expansion seemed to run off at a 45 degree angle (Figure below), which does not capture the true length of the polygon parcels.

An example output of what I am looking for is below.

This is a section where the original coastline line crossed the parcels and an "Intersect" was done. This new intersect polyline is symbolized by color (e.g. green = public). With this, the length of each segment can be totaled for the entire coast and used to determine percent ownership.

Any ideas on how to keep the raster expansion bounded to the parcel dimensions, or other ways about going at this will be extremely helpful! These images are a subset of one county - and I am looking to do this for 14 others. Using Arc 10.4.


texts: Mapnik - Preserve RGB value of transparent pixels

texts:
Currently generating PNG raster tiles (mapnik) with custom style where I store one information (map_feature) by channel (rgba).
Mapnik seems to multiply alpha with value of red, green or blue. I completly understand the rendering engine case seeing 0 as alpha interpreting as full transparent pixel. I just need it to write the color value anyway.
More explanations on my setup :

I have an OSM Database to deliver PNG raster tiles.
Nothing much to say about the setup - It's working perfectly.

PostgreSQL
Mod_tile
Mapnik

I've setup a new style based on osm-bright in order to display one type of information by channel :

Channel Red   rgba(255,0,0,1) --> Rivers
Channel Green rgba(0,255,0,1) --> Water_areas
Channel Blue rgba(0,0,255,1) --> Forest

Using alpha as 1 and with 3 informations works.
However when i try to swap alpha channel to 0 and add info in alpha channel  :

Channel Alpha rgba(0,0,0,255) --> Roads
Channel Red   rgba(255,0,0,0) --> Rivers
Channel Green rgba(0,255,0,0) --> Water_areas
Channel Blue rgba(0,0,255,0) --> Forest

When I'm analyzing the raster tile (the BINARY not with my eyes) I can't find anymore the value store where as rgb were.
I don't need to see them just disable the interpretation of Mapnik (or agg_renderer).
I've dig deep inside mapnik github code I found this :

 Github - Mapnik: test/unit/imaging/image_premultiply.cpp

It seems to be a function called set_premultiplied_alpha() which could explain why (255,0,0,0) pixel is interpreted as (0,0,0,0).

How can I fix my problem?


texts: Determining if polygon contains another smaller polygon using geometry objects in ArcPy?

texts:
I have a layer of parcels where some "parcel-itos" are inside larger "container" polygons, but these parcelitos are inside holes of the containers (e.g. they do not overlap, but they do touch). Both parcelitos and containers are in the same feature class. I would like to be able to relate the parcelitos' IDs to their containers' IDs. This is what I've tried:

inside = <Polygon object of Parcelito>
outside = <Polygon object of Container>

outside.contains(inside, 'BOUNDARY') #returns false
inside.within(outside, 'BOUNDARY') #returns false

From what I could glean on ESRI's info for Polygon Objects and Methods, it seems like specifying 'BOUNDARY' should have made both of the tests above return True. I don't want to use inside.touches(outside) because the iterative process that contains this test would also pick up the parcels surrounding the container polygon, and I don't care about those relationships.

Any idea of what's going on with these polygon object methods?


texts: Seeking algorithm behind End Point Rate (EPR) to calculate movement of shoreline over period of years?

texts:
I am analyzing the shoreline change over a period of 20 years(1998-2018). I have following data:

digitized shorelines of 1998, 2003,2008,2013 and 2018
one common baseline for all shorelines
transect lines from baseline to shorelines.
I have the lat/long of the intersection points of transect lines and
shorelines.
The distance between shoreline and baseline from the intersection
point.

I need to calculate the EPR (End point rate) and NSM (net shoreline movement).
My DSAS Tool is not working due to some .dll file missing.
I need to figure out the algorithm behind these calculations. Any suggestions?


texts: Only perform Near analysis on selected rows

texts:
I am writing a Python script to perform a near analysis on features from two different shapefiles with common attributes. This is happening in a Search Cursor.
However, whenever I tried to use the Near tool after I selected the relevant features from each shapefile OR created new layers of only the relevant features, the Near tool performs correctly on the selected rows but fills the rest of the table NEAR_FID and NEAR_DIST rows with -1. I have double and triple checked this process out of the cursor loop and my selection code works correctly. Obviously, this means that I cannot do this iteratively as each use of the tool will overwrite the previous use with -1 values.
I am currently solving this issue by using the Generate Near Table tool and then appending all tables together to join back to the original shapefile, but this takes a much longer time than the Near tool alone. According to the ArcGIS blog post Near By Group this tool should work in iterations. Is there something unique about the ModelBuilder iterator tools that allows this that doesn't exist in python Cursors?

Here is my cursor with the Near tool:

finishList = []
arcpy.MakeFeatureLayer_management("myGDB.gdb\\lb07", "log")
arcpy.MakeFeatureLayer_management("myGDB.gdb\\vms07", "vms")
rows = arcpy.SearchCursor("myGDB.gdb\\lb07")
for row in rows:
    id = row.fmc_logves_id
    date = str(row.string_date)
    finished = str(id) + " " + str(date)
    if finished in finishList:
        continue
    else:
        arcpy.SelectLayerByAttribute_management("log", "NEW_SELECTION", """"fmc_logves_id" = """ + str(id) + """ AND "string_date" = '""" + str(date) + """'""")
        arcpy.SelectLayerByAttribute_management("vms", "NEW_SELECTION", """"fmc_ves_id" = '""" + str(id) + """' AND "string_date" = '""" + str(date) + """' AND "activity" = 1""")
        arcpy.Near_analysis("vms", "log")
        finishList.append(finished)
del row, rows


texts: GDAL/OGR losing OSM closed way features that OGR2OGR imports

texts:
I have managed to import an OSM PBF file into QGIS using ogr2ogr -> SQLite and then importing layers using the GUI- (essentially as per the first step in this answer :

ogr2ogr -f "SQLite" -dsco SPATIALITE=YES texas.db texas-latest.osm.pbf

My version of GDAL and OGR2OGR is 2.4.1 (as reported by both - I checked)

That appears to work fine. I'm trying to duplicate this in a Python plugin without the SQLite intermediate step. Points, lines, and multilines appear to be readable fine (exact same numbers of features). However, the total number of multipolygon features is significantly less. Further investigation shows the missing features are things like buildings which have an osm_way_id set. I.e. these have been encoded in the PBF file as closed ways.

Here is an example imported with ogr2ogr:

The GDAL library code does not see any features with osm_way_id set! (although it recognizes this is a valid field name). Here is my code:

    gdal.SetConfigOption("OGR_INTERLEAVED_READING", "YES")

    input_pbfname = "C:/Data/QGIS/Texas/texas-latest.osm.pbf"
    ds = gdal.OpenEx(input_pbfname, gdal.OF_VECTOR)
    if ds is None:
        QMessageBox.information(None, "open failed", "Open failed")

    nly = ds.GetLayerCount()

    #version_num = int(gdal.VersionInfo('VERSION_NUM'))
    #QMessageBox.information(None, "Version",  str(version_num))

    nPoints = 0
    nLines = 0
    nMultiLines=0
    nPolygons = 0

    # GetNextFeature handles interleaved reading correctly
    for ily in range (0,nly):
        layer = ds.GetLayerByIndex(ily)

        layer.ResetReading()
        feat = layer.GetNextFeature()
        while feat is not None:
            if layer.GetName() == "points":
                nPoints = nPoints+1
            elif layer.GetName() == "lines":
                nLines = nLines+1
            elif layer.GetName() == "multilinestrings":
                nMultiLines = nMultiLines+1
            elif layer.GetName() == "multipolygons":
                nPolygons = nPolygons+1
                idx = feat.GetFieldIndex("osm_way_id")
                if idx>=0:
                    wayid = feat.GetFieldAsString(idx)
                    if (wayid is not None and wayid!="" and wayid!="0" ):
                        nWayId=nWayId+1

            idx = feat.GetFieldIndex("osm_way_id")
            if idx>=0:
                # tried both GetFieldAsInteger and GetFieldAsString
                # also tried counting all multipolygon features with osm_way_id set
                wayid = feat.GetFieldAsInteger(idx)
                if wayid == 646651935:
                    QMessageBox.information(None,"646651935 building found","Layer:"+layer.GetName() )

            feat = layer.GetNextFeature()

    QMessageBox.information(None, "Count", "Points: %d Lines:%d MultiLine:%d  Poly:%d" % (nPoints,nLines,nMultiLines, nPolygons) )

    # close the temporary layer
    ds = None

So why are these "closed ways" getting imported correctly by ogr2ogr but not the GDAL api? Is there a flag or setting that I'm missing? Or another subtlety with the interleaving that I have missed?  (I'm aware of the reasons for this with OSM data - in the past I've coded an OSM PBF reader in C#...)

ogr2ogr is importing 3385834 multipolygons, but the gdal API is only importing 10896. The other layers (point, line, multiline all match)


texts: OGR: "ERROR 1: No field definitions found for '...', is it a table "

texts:
I am writing a script to test if the number of elements is equal in the input file of a process and in the output (database). This is my code:

def check_output(value, refcount, reftype):
    db, idsch, table = value.split('.')
    ident, schema = idsch.split('_')

    dsn = ogr.Open(get_connstr())
    if dsn is None:
        raise Exception("Reading data failed.")
    lyr = dsn.GetLayerByName('{}.{}'.format(idsch, table))
    if lyr is None:
        raise Exception("Could not find the layer.")
    count = lyr.GetFeatureCount()
    if count != refcount:
        raise Exception("The number of elements written to the database is different from the number of elements in the input file.")

I get the following error:

ERROR 1: No field definitions found for 'buff_out', is it a table?
raise Exception("Could not find the layer.")
Exception: Could not find the layer.

The result of this (printing the argument of the GetLayerByName function):

print('{}.{}'.format(idsch, table)):

is this:

process-one-output_0489e558-eb28-11e7-a0d8-54ee758a7423.buff_out

If I understand it correctly, it sucessfully connects to the database but it cannot access the layer. Anyone knows what causes the error and how to solve it?


texts: How to delete a feature class from a GDB using ArcObjects for Java?

texts:
I've seen this solution which may works in ArcObjects for VB/.Net:

How to overwrite a shapefile or feature class using arcobjects

When I try to convert it into Java, I get a ClassCastException on the IDataset cast.

Here is my code:

IWorkspace workspace = ...
IWorkspace2 workspace2 = (IWorkspace2)workspace;
if(workspace2.isNameExists(esriDatasetType.esriDTFeatureClass, featureClassName)) {
    IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspace;
     // ClassCastException here
    IDataset featureClassDataset = (IDataset)featureWorkspace.openFeatureClass(featureClassName);
    if(featureClassDataset.canDelete()) {
        featureClassDataset.delete();
    } else {
        // alert: the feature class exists but couldn't be deleted
    }
}
IFeatureClass featureClass = featureDataset.createFeatureClass(...);

Does anyone know how we can't use the IDataset interface for a feature class as in VB?

How to delete a feature class in Java else?


texts:
I have a shapefile that looks like this:

I want to delete the holes, but I can't use any of the editing tools apparently...

What can I do?


texts: Finding areas of overlap using ArcGIS Desktop

texts:
I am going to have about 300 shapefiles representing the range of different species in the wild. What I want to be able to do is find areas where the most species occur to highlight areas where many of these species exist.

Is there a simple way to do that in ArcGIS Desktop 9?


texts:
is there a way to reproject a feature collection into a new one with another coordinate reference system? I am trying to achieve this by using geotools's reprojector:

SimpleFeatureSource featureSource = ...
SimpleFeatureCollection sfc = new ReprojectingFeatureCollection(featureSource.getFeatures(),
    org.geotools.referencing.CRS.decode("EPSG:4326"));

sfc presents another schema's CRS but its geometries looks the same before reprojecting. What am I missing? Do I have to manually create a new geometry with desired CRS and coordinates for each geometry in featureSource or geotools offer a more elegant and efficient way to do this? I am using geotools version 17.2.


texts:
I am trying to randomly distribute points in a raster (that has 11 unique values). I followed this posting and was able to recreate it, being partially successful.

As mentioned about, my raster has 11 unique values and I want different number of points randomly placed based on the unique strata. Also, I need different number of points.

So as the posting says,I isolated individual unique layers based on the unique value

raster[raster != 81] = NA # creating a raster with pixels of value 82 and remaining pixels Nodata

Then I used the sampleStratified() function as mentioned int he post

x<-sampleStratified(raster,1 , sp=T)

I want only 1 point to be randomly placed and I get this error

Error: size > 0 is not TRUE
In addition: Warning messages:
1: In .local(x, size, ...) :
fewer samples than requested found for stratum: 81

I don't understand how to rectify this because there are over 2000 pixels in the raster with pixels of value 81. So does this mean that sampleStratified() was not able to place 1 point in any of the >2000 pixels present? How do I rectify this error.


texts: How to create line segments between pairs of pixels in a raster?

texts:
I want to draw several lines between pixels (please see the image below) according to some criteria :

First Criterion: The distance between 2 pixels is less than 2000m and greater than 500m
Second Criterion: The absolute difference of Z values is less than 600m and greater than 400m

For example, let's consider two pixels of a DEM: I want to draw a line between their centroids if they satisfy both criteria.

I tested by converting my raster DEM to points and applying these criteria but it takes a lot of time.
Have you any proposal or suggestions?
I'm using GRASS, Postgis or Python


texts: Exporting geodatabase using R-Bridge with ArcGIS Pro

texts:
I am trying to export a spatial points data frame as a Arc geodatabase using arc.write function in the R-bridge (arcgisbinding) in R. I'm working in R v3.6.1, ArcGIS Pro v2.4.0 Advanced License, and arcgisbinding 1.0.1.241. The line of code will run for about 10-15 minutes before it fails and throws the following error:
> arc.write(path = "C:/Desktop/Geodatabase.gdb/fc", data = LLcoor, overwrite = T)
Error in .call_proxy("arc_write", path, pairlist(data = data, coords = coords,  :
  insert row failed

What's strange to me is that, about a week ago I could export these gdb's with no problems. The dataframe is quite large 350,000 rows and plan to add more. I checked for and removed all na's in my coordinates before creating the SPDF. Spatial Points Data Frame below:

I've tried un and reinstalling arcgisbinding, switching my write location, switching computers and clearing my internet cache.
Any thoughts?


texts: Splitting .tif image into several tiles?

texts:
I have an image of size 1GB (.tif), with the width and height 94000x71680. I would like to chunk this image into 20000X20000 tiles so that I can process them.

How can I do this?


texts: How to transfer data from spreadsheet into shapefile?

texts:
I am looking at Deprivation in the UK and there is some good data in the form of spreadsheets on this website http://www.communities.gov.uk/publications/corporate/statistics/indices2010

I would like to know how i can transfer some of the data into a shapefile??


texts: We need to link to Quantum GIS on a major official website. My legal department asks me if "Quantum GIS" is trademarked. Is it?


texts:
I have a table of tessellated hexagons that I created as a hex binning grid to cover an area.

I've then overlaid a point set to give some of the hexagons a value, see diagram:-

The shapes with value "0" have no value and I'd like to interpolate these from the values of the neighbours as per Kriging. How can I accomplish this in postgres given my table looks like:-

hex_grid_data
(
 gid SERIAL,
 wkb_geometry (geometry),
 value INTEGER
)


texts: "AttributeError: 'NoneType' object has no attribute 'fetchall'" loading data to Elastisearch

texts:
I want to load data from PostgreSQL to Elasticsearch and I write below code, but it gives me the error in for loop, my code is as below:

from elasticsearch import Elasticsearch

#res = requests.get('http://localhost:9200')
es = Elasticsearch([{'host': 'localhost', 'port': 9200}])
sql_conn = psycopg2.connect(host="localhost",database="priyanka",user="postgres",password="admin")
sql_cursor = sql_conn.cursor()
sql=('SELECT * FROM "Sales"')
sql_cursor=sql_cursor.execute(sql)
tempList=[]
tmp = {}
for row in sql_cursor.fetchall():
    print("1")
    Sales_id =str(row[0])
    Product_name =str(row[1])
    sales =str(row[2])

    tmp["Sales_id"]=Sales_id
    tmp["name"]=Product_name
    tmp["sales"]=sales

Error is:

Traceback (most recent call last):
  File "D:/Priyanka Pagar/Softwares/Team/Python37/elastic - Copy.py", line 16, in <module>
    for row in sql_cursor.fetchall():
AttributeError: 'NoneType' object has no attribute 'fetchall'


texts: Decreasing Rows and Column of Raster file?

texts:
I am working on urbanization modelling.

While processing i need to have the columns and rows of the raster file less than 100. But my raster file has more than 2000 rows and columns. So, I re sampled it but after that the clarity was too poor.

What can be the other techniques to reduce its number without compromising too much of clarity?


texts: Finding Gaps in Hundred Range Street Centerline Segments in GIS

texts:
Is there a way to find gaps for Hundred Range Street Centerline Segments within the attribute table?

An example that I can think of, is if there is a hundred range street segment that shows 1001-1049 W Adams St and the following centerline segment shows 1071-1099 W Adams St. The missing data can effect our composite address locator and create other issues when there is a gap in the hundred ranges. I

f there was a capability to find which centerline hundred ranges are not continuous within our database it would be useful.


texts: Heatmaps: right terms (isopleth vs choropleth, etc) and how to implement a value-based one that looks good?

texts:
NOTE: My original question was about interpolation methods, but I've found that I was actually asking the wrong question. Interpolation methods with Heatmaps and IDW are the least of my problems right now.

I've been investing a lot of time trying to understand heatmaps visualizations in the context of GIS in the last month. There's a lot of information, theory and a gazillion different implementations and stacks, it's very easy to be lost in this sea of information, hence, I'm here asking for help ;)

Right now, I'm using PostGIS + GeoServer as map-sever middleware + Leaflet for the client.

But my main question is really about the visualization below:

https://www.walkscore.com/apartments/search/CA/Santa_Monica

Screenshot:

And also, here's one of the tiles:

1) Is this a choropleth or a heatmap, or what is the right term for this visualization? If someone knows specifically what kind of heatmap that is, please let me know.

2) I'd love to emulate this style with geosever (or any other software / stack of software for tha matter, this one is served as tiles from the server). Should I just invest more time learning SLD or should I just go custom and implement the rendering code myself? Does anyone with more experience than I have put in me in the right direction?

The data

Lastly but not least, the data that I'm plotting is simple, but it might be tens of thousands of points representing buildings/apartments/houses in a city, and the heatmap would be value-based, by revenue per month, for instance. So, lat/lon + a value (revenue). I want to plot that in a heatmap similar to the one shown above.


texts:
What exactly is NoData? Since you can set NoData. Is it just a rendering after reading the data and finding areas where there literally is no data but there should be? For example a raster is constrained to 500x500, but the data only exists for 443x300, so during the read, it fills out the remaining area without data with whatever NoData is set as? Meaning QGIS would fill it out differently than ArcMap which would fill it out differently than GDAL since it happens during read? Or is NoData actually hard saved as those numbers in the file?
How can I access no data in GDAL? On the raster I'm working on, it seems to be in e notation 3.402823...e+38.  If I wanted to do something like below with the band array to check for no data values, what is the appropriate thing to put in the if statement?
band = ds.GetRasterBand(1)
band_one = band.ReadAsArray()

for row in band_one:
    for column in row:
        if column is [*no_data*]:
            print("No data value found")


texts:
I've been searching for a good example of how to publish a Web Process Service (WPS) on the GeoServer platform.  The page where you'd think they'd document this process is actually pretty thin.

There's some documentation on implementing an OWS service, but I am unsure what the conceptual and implementation overlap is between OWS and WPS in GeoServer.

Is anyone aware of a good example of implementing a WPS on GeoServer?


texts: Can't join attributes by location between 2 polygon layers in QGIS 2.18

texts:
I have a shapefile of Spanish municipalities downloaded here and another shapefile of functional urban areas (FUA) downloaded here.

I want to add to each municipality the urban area id if the municipality is inside a urban area. To do so, I go to the menu Vector > Data management tools > Join attributes by location with the following parameters :

target layer: Municipios
join vector layer: FUA
predicate: all checked
take attributes of the 1st located feature selected
Keep all records selected
Create a temporary layer selected

The columns of the joined layer which should contain the values from the FUA layer are all empty.

I've also checked the geometry on both layers, the municipalities one has 4 invalid polygons and even removing them, it is still not working.

What am I doing wrong ?

I had no problem before joining attributes on QGIS 2.14.


texts: How do I use custom fonts for labelling in Geoserver?

texts:
I want to use a custom font I have downloaded from the web with Geoserver for my labels. How do I specify the location of my font so that the font-family css parameter can read it?

Geoserver 2.1.4


texts: Return table names for query layers using arcpy

texts:
The dataSource property for layers does not return the name of the tables used in query layers. Instead, it returns something like the following:

'Database Connections\\ConnectionFile.sde\\sde.SDE.%queryLayer'

Whereas the tables in SDE that the query layer is referencing may be called something like sde.SDE.FEATURECLASS (and this is what dataSource would return if the layer were a non-query layer layer).

Since you can't use the dataSource property to return the tables used to make a query layer, how do you return the table names?

Using ArcGIS/ArcPy 10.3 and ArcSDE 10.1


texts: Align raster to fishnet using ArcGIS 10.3

texts:
I am having some trouble aligning a raster (elevation) to a fishnet I created. The cell sizes for both the raster and the fishnet are 30x30 meters. I have tried using the clip and extract by mask tools to clip the raster to the fishnet; I made sure that the extent was set to match that of the fishnet. Any thoughts?

For those asking for a screenshot for added clarity. The first image is an example of another raster that I was able to successfully align to the fishnet. In this case, I converted it from a polygon to a raster and made sure the output coordinates and processing extent were the same as my fishnet.

The second image is my elevation raster that will not align to the fishnet.


texts: QgsMapCanvas window flashes on then disappears

texts:
I am trying to make a standalone app, to be able to select and download some spatial data that I keep around.
I've read around and found I can do that with Qgis and Qt, so I grabbed some example code from the docs, and assigned it a shapefile to open.
And it wouldn't work, regardless of what I tried with PATH and PROJ_LIB.
While I was trying to work out my issues with:
proj_create_from_database: Cannot find proj.db
proj_get_authorities_from_database: Cannot find proj.db

I finally found out, I can avoid that by using:
C:\Program Files\QGIS 3.16\bin>python-qgis-ltr.bat myApp.py

Now a window is shown for a split second, then it disappears and the program finishes, as cmd is ready for new inputs.
Here is the whole code:
import sys
sys.path.append(r"C:\Program Files\QGIS 3.16\apps\qgis-ltr\python")
import os

os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = r'C:\Program Files\QGIS 3.16\apps\Qt5\plugins'
os.environ['PROJ_LIB '] = r'C:\Program Files\QGIS 3.16\share\proj'
os.environ['PYTHONPATH'] = r'C:\Program Files\QGIS 3.16\apps\Python37'

from qgis.PyQt.QtGui import (
    QColor,
)

from qgis.PyQt.QtCore import Qt, QRectF

from qgis.core import (
    QgsVectorLayer,
    QgsPoint,
    QgsPointXY,
    QgsProject,
    QgsGeometry,
    QgsMapRendererJob,
)

from qgis.core import *

from qgis.gui import (
    QgsMapCanvas,
    QgsVertexMarker,
    QgsMapCanvasItem,
    QgsRubberBand,
    QgsMapToolPan,
    QgsMapToolZoom,
)

from qgis.PyQt.QtWidgets import QAction, QMainWindow

class MyWnd(QMainWindow):
    def __init__(self, layer):
        QMainWindow.__init__(self)

        self.canvas = QgsMapCanvas()
        self.canvas.setCanvasColor(Qt.white)

        self.canvas.setExtent(layer.extent())
        self.canvas.setLayers([layer])

        self.setCentralWidget(self.canvas)

        self.actionZoomIn = QAction("Zoom in", self)
        self.actionZoomOut = QAction("Zoom out", self)
        self.actionPan = QAction("Pan", self)

        self.actionZoomIn.setCheckable(True)
        self.actionZoomOut.setCheckable(True)
        self.actionPan.setCheckable(True)

        self.actionZoomIn.triggered.connect(self.zoomIn)
        self.actionZoomOut.triggered.connect(self.zoomOut)
        self.actionPan.triggered.connect(self.pan)

        self.toolbar = self.addToolBar("Canvas actions")
        self.toolbar.addAction(self.actionZoomIn)
        self.toolbar.addAction(self.actionZoomOut)
        self.toolbar.addAction(self.actionPan)

        # create the map tools
        self.toolPan = QgsMapToolPan(self.canvas)
        self.toolPan.setAction(self.actionPan)
        self.toolZoomIn = QgsMapToolZoom(self.canvas, False) # false = in
        self.toolZoomIn.setAction(self.actionZoomIn)
        self.toolZoomOut = QgsMapToolZoom(self.canvas, True) # true = out
        self.toolZoomOut.setAction(self.actionZoomOut)

        self.pan()

    def zoomIn(self):
        self.canvas.setMapTool(self.toolZoomIn)

    def zoomOut(self):
        self.canvas.setMapTool(self.toolZoomOut)

    def pan(self):
        self.canvas.setMapTool(self.toolPan)

qgs = QgsApplication([], True)

qgs.initQgis()

pathToShp = r"D:\Docs\myTestShape.shp"
inLayer = QgsVectorLayer(pathToShp, "myTestShape", "or")

w = MyWnd(inLayer)
w.show()

What do I need to do to see something displayed and pan around it?


texts:
I am currently trying to choose suitable fonts for a map. A very time consuming experience - actually too time consuming as i´m not a trained typographer and will inevitably make all the classic mistakes!

In order to save myself time, and probably many other GIS users who very probably make all the classic errors when choosing font types, does anyone have any SIMPLE guidelines regarding which type of fonts should be used for standard thematic maps. I am not talking about artistic fonts for special types of graphic, just simple clear easy to read fonts which work well together in different sizes and weights in the same map.

My main problem at the moment is trying to choose a large bold clear font for the main point theme layer annotations, and a less-weighted font for the place names. I want to have a sans-serif for the point layer which displays the names of the points in a dominant concisive factual way, but want a fitting and complementary serif-font for the place names which is thinner and less dominant.

If anyone has standard fonts which they use together for most of their work, I would be grateful for any advice,


texts:
I am wondering what the purpose of arcpy.Array() is. It seems an intermediate step between creating a list of objects and performing some sort of conversion, for example, such as in the sample code for the 10.1 help.

    features.append(
            arcpy.Polyline(
                arcpy.Array([arcpy.Point(*coords) for coords in feature])))

What purpose does the arcpy array serve? Why can't a list suffice? What's the logic behind this intermediate step?


texts: Install ArcGIS Desktop unattended or silently

texts:
I’m using the instructions for Installing ArcGIS Desktop silently. However the process doesn’t always work in unattended mode, because of a prompt about "The following applications need to be closed before continuing installation: ArcGIS Cache Manager". What switch or parameter is needed to auto-select an answer and keep going?

The current command line:

%SystemRoot%\System32\msiexec.exe /I ^
    "%cd%\ArcGIS_Desktop\Desktop\SetupFiles\setup.msi" ^
    ESRI_LICENSE_HOST=LICENSESERVER ^
    SOFTWARE_CLASS=Viewer ^
    SEAT_PREFERENCE=Float ^
    INSTALLDIR=C:\ArcGIS ^
    INSTALLDIR1=C:\Python27 ^
    DESKTOP_CONFIG=TRUE ^
    ADDLOCAL=ArcMap,ArcCatalog,Python,SpatialAnalyst ^
    /qb


texts: Can't split polyline into 20 metre segments

texts:
I'm trying to make 20 metres segments out of one line, but I'm getting this error:

  segment outside line, no segment created

I'm using the settings shown in the screenshot. If I unblock the "force segments to be exactly of given lenght(...)" it works but then I get segments of 30 metres instead. The segments need to be 20 metres.


texts: Implementation for Boolean operations on polygons like union, intersections, contains?

texts:
https://en.wikipedia.org/wiki/Boolean_operations_on_polygons#Algorithms
I wish to understand how are Boolean operations on polygons implemented in GIS related geometry libraries. I've tried to go deeper into Shapely, OGR, JTS in that chronology. I believe that the former is a port of the latter, in all the three libraries. I advanced until here, and got stuck. It seemed to me that this function is called recursively (I'm sure I'm wrong here).
Could anyone here who already has contributed to any of the libraries explain?


texts: OpenLayers-3 setting zoom so layers fit in map window

texts:
Is there a way to dynamically set the zoom level so that the current layers always fit within the width of the map frame on page load?

Currently the zoom level is set statically such that when using a full resolution monitor the application loads perfectly the entirety of the map fits but when loading the map at a lower resolution the map loads at the same zoom level but features are cut off.


texts: Select "interior" lines of a polygon-to-line layer

texts:
I have a polygon featureclass that has adjacent features. Using Arc10, I am turning those polygons into lines using the "Featureclass to Lines" tool.

Now I would like to select ONLY THE INTERIOR line(s) (color=grey) so that I can run some ADDITIONAL processes on them. The key is that some (but not all) feature may touch the boundary line (color=blue). I am trying to use "Select By Location", but I cannot come up with a way to accomplish this task.

I might just be overthinking this one, but I am stuck. Any suggestions?


texts: Is it poor form to rotate Albers US state projections if you're only viewing a single state?

texts:
We are using the Albers US map in a project with Alaska and Hawaii shown as well. The project allows users to zoom into a state and only see that state for more detail related to the data we have. When zoomed into the state, we felt the state displayed better when rotated to view the single state in a more straight on view. For example Oregon and Washington would rotate counter clockwise so the top and bottom are aligned with horizontal view.

Is this fine to do? And how is the same rotation viewed if done within the US projection as well when zoomed into a state that you can still view surrounding states?


texts: TopologyException error during fortify SpatialPolygonsDataFrame

texts:
I tried to reproduce a recent question on spatial data visualization using ggplot2 and received the following error message during fortify(world, region = "ISO3"):

Error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point
3980740.7660219339 4277735.2740812832 at 3980740.7660219339 4277735.2740812832

Any ideas on that?


texts: QGIS - create shapemap with region polygons from flatfile containing us zip to region assignments

texts:
I'm new to QGIS and I hope that you can help me with my issue. I have a csv file containing the US zip codes, their latitude / longitude and a Region (Q1 - Q6). For my Dashboard tool I'd like to create a shapefile in GeoJSON, where the zip codes are merged into the 6 regions.

Is there any easy way to achieve this in QGIS, e.g. create layers based on a specific attribute?


texts: Choosing correct local coordinate system?

texts:
I need to use GIS data from various parts of the world. As a result GIS data is in different coordinate systems.

How to find out the correct local coordinate system?

Is there the official list of them?

For example, I should choose local coordinate system for Chicago, Ibadan (Nigeria) and Rio De Janeiro.


texts: Efficiently classify raster with irregular polygons

texts:
I have a project on my desk to classify an elevation model using different elevation breaks for each polygon in the area of interest.

In pseudo code:

FOR each polygon IN layer:
    CREATE_MASK FROM this polygon
    CLASSIFY elevmod WITH this polygon's LOOKUP_TABLE
    SAVE_AS elevmod_classified_this_polygonID

MERGE elevmod_classified*
CONVERT from 16bit float to 256 index colour palette

What is the most efficient way to do this?

The workflow I have now is wasteful of processing time and storage. Each classify step creates a new raster file for every polygon, and then all are combined into one at the end.  The source raster is 650mb compressed and about 30 zone polygons (so 30 times some portion of 650mb for extra storage). This seems rather silly since I can classify a whole raster at run-time in a second using the standard symbology tools.

Plus the elevation breaks for each polygon haven't been finalized, so I'll need to run through this exercise a number of times before everyone is satisfied.

The classification elevation lookup tables resemble the following. Each number is the upper elevation boundary of that class:

polygonID, class1, class2, class3, ...
169,  850, 1100, 1300, 1450
170, 1000, 1250, 1400, 1600

The final classified raster has cell values of 1,2,3,4,5... corresponding to class1, class2, ... Though this too is open for change.

I have Qgis 1.7 and Arcinfo 10 with Spatial Analyst at my disposal (and Arcview3 for that matter).


texts: Load partial shapefile into Postgis / GeoDjango project

texts:
I have a shapefile with Canadian postal codes, but am only looking to load a small subset of the data. I can load the entire data file and use SQL or Django queries to prune the data, but the load process takes about 2 hours on the slower machines I'm using.

As the data I'm actually after is about 10% of the dataset, this isn't a very efficient process.

I'm following the instructions in the Geodjango tutorial, specifically the following code:

from django.contrib.gis.utils import LayerMapping
from geoapp.models import TestGeo
mapping = {'name' : 'str', # The 'name' model field maps to the 'str' layer field.
               'poly' : 'POLYGON', # For geometry fields use OGC name.
               } # The mapping is a dictionary
lm = LayerMapping(TestGeo, 'test_poly.shp', mapping)
lm.save(verbose=True) # Save the layermap, imports the data.

Is there a way to only only import data with a particular name, as in the example above?

I'm limited to the Linux / OS X command line, so wouldn't be able to utilize any GUI tools.


texts:
This question was stimulated by a discussion between WhiteBox and myself.

Both Reclass and Con (i.e. a conditional evaluation tool) are examples of tools that implement a Switch statement or nested If-Then-Else statements. But when should one use Reclass and when is it more appropriate to use Con (i.e. a conditional evaluation tool)?

I generally think of reclass as a simple tool and Con as a more complex tool and whiteboxdev the other way around.

I use Con mostly when I want to pull data from another raster layer based on the condition. I use reclasss when all the info is in the origin raster. I use con to "pick" values to replace in a Raster and reclass to do a "suite of changes"

Any best practices or what do you do?


texts: ECEF (EPSG:4978) glTF model doesn't show up in Cesium

texts:
The issue:

I have a glTF (ver. 2.0) mesh (converted from ply to dae to glTF) with triangles coordinates in EPSG:4978 that doesn't show up in Cesium (ver. 1.42).

The code:

This is the code in the main js app:

var model = scene.primitives.add(Cesium.Model.fromGltf({
    url : 'meshes/MESH_4978.gltf', // MESH_3857.gltf
    show: true,
    modelMatrix : Cesium.Matrix4.IDENTITY,
    scale : 1.0,
    debugWireframe: true
}));

The glTF file itself is well rendered here: https://gltf-viewer.donmccurdy.com/ (but the horizontal plane is vertical at start).

And I can see that the glTF file is well loaded with a Status=200 in the Firefox browser network monitor (or in Chrome).

The question:

So, where am I wrong?
Which coordinates system does Cesium use?

Other try:

I also tried with mesh coordinates in EPSG:3857; the exact same thing happens, i.e. file is loaded but nothing shows up.


texts: Can MapInfo Professional use data direct from a File Geodatabase?

texts:
Esri released its File Geodatabase API in mid 2011, and it seems that QGIS via GDAL's OGR FileGDB Driver added support for the File Geodatabase within about a year.

By support, I mean simply being able to display, identify and query features using a direct read, with editing them not necessarily being required.

Has such support of File Geodatabases also been implemented for MapInfo Professional at any release (or via an extension)?

Please note that this question is only about File Geodatabase (*.gdb) rather than Personal Geodatabase (*.mdb) support because there is a big difference between the two formats.


texts: Converting CIB data to Google Tiles with GDAL

texts:
I'm trying to convert CIB map data into Google compatible Tiles through the following process (in Linux, GDAL compiled from source):

gdalbuildvrt data.vrt A.TOC
gdal_translate -of VRT -expand rgba data.vrt data2.vrt
mkdir data
gdal2tiles.py data2.vrt data

gdal2tiles.py gives the following error:

  Generating Base Tiles: ERROR 1:
  VRTSourcedRasterBand::IRasterIO()
  called recursively on the same band.
  It looks like the VRT is refe
  rencing itself. ERROR 1: IReadBlock
  failed at X offset 0, Y offset 0 ERROR
  1: VRTSourcedRasterBand::IRasterIO()
  called recursively on the same band.
  It looks like the VRT is refe
  rencing itself. ERROR 1: IReadBlock
  failed at X offset 0, Y offset 0
  Traceback (most recent call last):
  File "/usr/local/bin/gdal2tiles.py",
  line 2241, in ?
      gdal2tiles.process()   File "/usr/local/bin/gdal2tiles.py", line
  478, in process
      self.generate_base_tiles()   File "/usr/local/bin/gdal2tiles.py", line
  1276, in generate_base_tiles
      dsquery.WriteRaster(wx, wy, wxsize, wysize, data,
  band_list=list(range(1,self.dataBandsCount+1)))   File
  "/usr/lib/python2.4/site-packages/GDAL-1.8.0-py2.4-linux-i686.egg/osgeo/gdal.py",
  line 746, in WriteRa
  ster
      buf_pixel_space, buf_line_space, buf_band_space ) TypeError: not a
  string

I know my CIB data source is valid, am I missing an obvious step?  Is this a user error or a bug in the GDAL software?  Thanks so much for any help/opinions.


texts: How to use the math tool ("times") in ArcGIS modelbuilder to multiply two raster sets while using an iterator?

texts:
I have two sets of rasters (in the same gdb), representing "loss" and "forestmask".

I would like to multiply each forestmask raster (value 1 = forest) with its corresponding loss raster (value 1=loss) to create a raster that shows deforestation (in which value 1 = deforestation).

Each loss raster should be linked to one corresponding forestmask raster (they have the same extent):
e.g.

GFC2014_forestmask_00N_020E * GFC2014_loss_00N_020E = GFC2014_deforestation_00N_020E (=output filename)

I guess this should not be too hard to do, and I tried to do this with the ArcGIS modelbuilder, iterating through the "forest mask" rasters using a wildcard, but I got stuck when I wanted to make use of the Parse Path tool and in-line variable tricks for defining the second input file (i.e. the loss raster). Right now, I tried to add a new variable (=file) in the model, called loss below, with input:

"GFC_2014_loss_"+ "%Name%"[21:]

But this is not working: the Times tool does not accept the loss variable as it is as the second input. (not sure if the syntax is correct like this, and if this can be used in a "file" variable).

I think I know how to do this in R with the raster package, but I heard Python is faster for these kind of calculations, so any help is welcome. Preferably I would like to do this in the model builder so that I can easily link it to other models I created.


texts: Using zoo-project WPS service in QGIS WPS-Client

texts:
I've created a python wps service in the zoo-project framework, which is used to calculate a Line of Sight based in on a start point, an end point and a start point height. The service works fine if i send a URL like the following:

127.0.0.1/cgi-bin/zoo_loader.cgi?Service=WPS&Request=Execute&Identifier=ExtractlineofsightGML&DataInputs=startpoint=<gml:Point><gml:coordinates>711330,6197923</gml:coordinates></gml:Point>;endpoint=<gml:Point><gml:coordinates>711370,6197970</gml:coordinates></gml:Point>;height=20&Version=1.0.0&

I want to get the wps to work with qgis, but after connecting using the WPS-Client plug-in I first got an error saying that the inputs had to be gml (it was originally geojson). After changing my datatype to gml, as in:

 mimeType = text/xml
 encoding = UTF-8
 schema = http://schemas.opengis.net/gml/3.1.0/base/feature.xsd

it still doen't work. I can't choose a GML layer, which is in the layers panel as input, but I can chose a newly created scratch layer.

Unfurtunately when choosing the scratch layer i get the following error:

In my function I use psycopg2 to connect to a db and input a point as gml into the ST_GeomFromGML function

Anyone have an idea to where I'm failing?


texts: Adding multiple ground overlays to Google Maps API to toggle with checkmarks

texts:
I found this answer here by user @hoogw very useful regarding Overlaying ArcGIS REST API on Google Maps.
I've been trying to add several layers to be toggled by checkmarks to no avail. The base URL can be the same, but the layers will be different.
I know it's probably best to list all the overlays and push them to an array, then loop over them, but I'm too much of a novice to get it to work.
Below is the code referenced:
var map;
var Overlays= []
var Overlay;

function get_overlay_image(_ne, _sw){

                var imageBounds = {
                  north: _ne.lat(),
                  south: _sw.lat(),
                  east:  _ne.lng(),
                  west:  _sw.lng()
                };

            var export_bbox = _sw.lng() + ',' + _sw.lat() + ',' + _ne.lng() + ',' + _ne.lat();

        var url_param = '&format=png&transparent=true&f=image&bboxSR=4326&imageSR=3857';

         var _width = map.getDiv().offsetWidth;
         var _height = map.getDiv().offsetHeight;

         var _size = '&size='+ _width + ',' + _height

         var _layers ='&layers=show:4,5'

         var root_url =
     'https://ndgishub.nd.gov/ArcGIS/rest/services/All_GovtBoundaries/MapServer/export?bbox=';

     var image_url = root_url + export_bbox + url_param + _size + _layers;

    // delete remove all overlay in overlays array.
       while(Overlays[0])
                {
                Overlays.pop().setMap(null);
                }

     var  Overlay = new google.maps.GroundOverlay(
            image_url,
            imageBounds);

       Overlay.setOpacity(0.5)

       Overlay.setMap(map);

       // push new overlay into overlays array
       Overlays.push(Overlay);
}

function initMap() {

map = new google.maps.Map(document.getElementById('map'), {
    zoom: 9,
    mapTypeId: google.maps.MapTypeId.ROADMAP,
   // center: {lat: 38.907292, lng: -77.036420} //washington
    center: {lat: 48.16, lng: -103.61}
    //48.16, lng: -103.61
  });

google.maps.event.addListener(map, 'idle', function() {
              var bounds =  map.getBounds();
              var ne = bounds.getNorthEast();
              var sw = bounds.getSouthWest();

            get_overlay_image(ne, sw);
     });
  }

Maybe something like
overlayMaps = [
    { //layer1
    url:  {
    return url1 + layers ;
    },
    { //layer2
    url:  {
    return url2 + layers ;
    }]

function toggleLayer(i) {
  if (overlayMaps[i].getMap() === null) {
    overlayMaps[i].setMap(map);
  }
  else {
    overlayMaps[i].setMap(null);
  }
}

But how to implement with the various parameters?


texts:
I have a problem today with the Nominatim reverse geocoder for Leaflet.

I found a solution on github : https://github.com/perliedman/leaflet-control-geocoder , but I don't know how to use it.

I can use the mapzen solution for the normal geocoding but for the reverse it doesn't work...

In fact, I have a HTML page which contains a leaflet map at the right, and a toggle menu at the left. In this menu, I have a mathematic tool which make some treatments to find some Lat Long coordinates.
Those coordinates appears in a  with the dd.dddddddd format.

And I want this LatLong result to be located on my map with a point and a toggle popup.

I've made a little Schema to illustrate my needs :

I've tried to put a "L.popup" solution but it doesn't work :

L.popup();

function onMapClick(e) {
    popup
        .setLatLng(['sch-lat','sch-long'])
        .setContent("The coordinates are" + ['sch-lvpk-lat','sch-lvpk-long'].toString())
        .openOn(map);
}

map.on('click', onMapClick);

The 'sch-lat','sch-long' are the html inputs in which we have lat long values.

And I also tried to make a Nominatim reverse geocoding by using leaflet-control-geocoder but I don't know how to use it...
I made this :

L.control.showAll({bounds: L.latLngBounds(L.latLng(41.058, 9.257), L.latLng(51.009, -5.295))}).addTo(map);
L.control.geocoder('mapzen-vv6kj8z').addTo(map);
L.control.geocoder('Nominatim').addTo(map);
L.control.scale({options: {position: 'bottomleft',maxWidth: 100,metric: true,imperial: false,updateWhenIdle: false}}).addTo(map);
L.control.mousePosition().addTo(map);

But it doesn't work neither...

If you have another solution I accept it :D


texts: What is the right methodology to assign costs to edges of a network?

texts:
I am producing a routing facility for cycling safety based on an OSM network. The edges of this network have costs and these are calculated from different criteria both related to intrinsic (as the grade of road) and extrinsic (number of accidents per edge) characteristics. I already have found a method to assign costs but I do not know how to scale and weight them in order to make the routing coherent with my objective. My objective is to avoid dangerous edges assigning them higher costs and prefer safe roads and cycling paths, Dijkstras pg_routing actually does that but dangerous points are avoided "too much" and least cost paths are safe but incredibly long and "wriggly" . I tried to search some papers but I have not found anything related to my case, they are all much related to the criteria definition rather than to their application to a network model.

If someone could help I really would be very thankful.


texts: Storing datetime values in Shapefiles using QGIS?

texts:
I have a shapefile with date in a string field. I create a new field with date type and use the calculator to make the conversion. But the result only provide the date part, the time part is removed. Even when the preview show the correct conversion.

Here you can see the original string field time and the result date time2. And the use of the conversion function to_datetime() and the preview result at the end

If you go to the layer properties you can change the edit mode to date/time. Then if go to the data and try to edit time2 you get the current date/time. So is possible save date and time in the same field. But the update tool still doesn't work as expected

As J.R said. Even when you can change QGIS edit to allow Date/Time the save to shapefile will discard the time.


texts:
I am using arcpy to merge multiple feature classes from a list of feature classes. They all have the same fields. Let's say they all have the fields = ['CO_FIPS', 'MILES', 'FLOOD_ZONE', 'VALID', 'REGION']. However, I only want to include 'CO_FIPS', 'MILES',  and 'FLOOD_ZONE' in the output merged feature class. It seems I need to use field mapping to limit the fields in my output, but I am not 100% how to do that:

working_dir = r"C:\Projects\MyProj\stream_lines.gdb"
arcpy.env.workspace = working_dir
s_studies_for_merge = ['stream1', 'stream2', 'stream3', ... 'stream10']
fieldMappings = arcpy.FieldMappings()
merged_output = os.path.join(working_dir, "AllStreams")
fieldMappings.addTable('stream1')

for field in fieldMappings.fields:
  if field.name not in ['CO_FIPS', 'MILES', 'FLOOD_ZONE']:
   fieldMappings.removeFieldMap(fieldMappings.findFieldMapIndex(field.name))
arcpy.Merge_management(s_studies_for_merge, merged_output, fieldMappings)

My confusion is mainly that most examples show only 2 features being merged and that the mapping involves settings that I don't need for certain fields. I only want to exclude the fields I don't want and keep the fields I do. What do I need to do here?

UPDATE
I tried applying the below answer(@geodranic) this way (assuming I needed to iterate and add field mapping for each of the 10 feature classes in my s_studies_for_merge list):

mappings = arcpy.FieldMappings()
fields = ['REACH_ID', 'STUDY_ID', 'CID', 'FLD_ZONE', 'VALIDATION_STATUS', 'STATUS_TYP', 'MILES', 'STATUS_DATE', 'STUDY_TYPE', 'LINE_TYPE', 'BS_ZONE', 'BS_STDYTYP', 'PRELIM_DATE']
for item in s_studies_for_merge:
    print(item)
    for field in fields:
        print(field)
        map = arcpy.FieldMap()
        map.addInputField(item, field)
arcpy.Merge_management(s_studies_for_merge, s_studs, mappings)

I got an error that the merge failed because of one of the fields I didn't want included in the merge:

ExecuteError: ERROR 001156: Failed on input OID 1, could not write value 'VIII' to output field FEMA_Region
Failed to execute (Merge)

How can I apply this correctly?


texts: Updating raster color ramp and exporting new raster image using PyQGIS?

texts:
I have a .tiff raster layer and I want to update raster color ramp with new min/max values using PyQGIS and export to .png image with new color ramp.

The update raster color ramp I can do with code from this post :

from PyQt4.QtCore import *
from PyQt4.QtGui import *

layer = iface.activeLayer()

renderer = layer.renderer()
provider = layer.dataProvider()
extent = layer.extent()

ver = provider.hasStatistics(1, QgsRasterBandStats.All)

stats = provider.bandStatistics(1, QgsRasterBandStats.All,extent, 0)

if ver is not False:
    print "minimumValue = ", stats.minimumValue

    print "maximumValue = ", stats.maximumValue

if (stats.minimumValue < 0):
    min = 0

else:
    min= stats.minimumValue

max = stats.maximumValue
range = max - min
add = range//2
interval = min + add

colDic = {'red':'#ff0000', 'yellow':'#ffff00','blue':'#0000ff'}

valueList =[min, interval, max]

lst = [ QgsColorRampShader.ColorRampItem(valueList[0], QColor(colDic['red'])),
        QgsColorRampShader.ColorRampItem(valueList[1], QColor(colDic['yellow'])),
        QgsColorRampShader.ColorRampItem(valueList[2], QColor(colDic['blue']))]

myRasterShader = QgsRasterShader()
myColorRamp = QgsColorRampShader()

myColorRamp.setColorRampItemList(lst)
myColorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
myRasterShader.setRasterShaderFunction(myColorRamp)

myPseudoRenderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(),
                                                    layer.type(),
                                                    myRasterShader)

layer.setRenderer(myPseudoRenderer)

layer.triggerRepaint()

but how to export this raster to new .png image with new colors ?


texts: Converting dwg to shapefile using arcpy and ArcGIS js

texts:
The following code in arcpy is used to convert a dwg file into shapefile and it works perfectly when you run it within ArcMap:

import sys, zipfile, arcpy, os, traceback
from arcpy import env
# Set environment settings
env.workspace = "C:/GISApplication/GIS/Layer_Updated_Drawings"

# Set local variables
fileName = arcpy.GetParameterAsText(0)
featureType = arcpy.GetParameterAsText(1)
file = os.path.join(fileName, featureType)
inFeatures = [file]
outLocation = "C:/GISApplication/GIS/Layer_Updated_Drawings/Converted"

# Execute FeatureClassToGeodatabase
arcpy.FeatureClassToShapefile_conversion(inFeatures, outLocation)

However when I try to run this script using the following JavaScript code:

url = ARCGIS_SERVER + 'Toolbox/DwgToZippedShapefile/GPServer/LayerUpdate';
convParams = {
FileName: "MIPBoundary.dwg",  //Hardcoded for testing
FeatureType: "Polyline" //Hardcoded for testing
};
convExportToolName = "SHPFileResults";
convJobFolder = "DwgToZippedShapefile_GPServer";
var gpConversionProcess = new Geoprocessor(url);
gpConversionProcess.submitJob(convParams, gpConversionComplete, gpConversionStatus, gpConversionFailed);

it keeps failing and I keep getting the following error on arcgis server manager:

  Error executing tool. LayerUpdate Job ID:
  j7af4b50ca616491eb13dce08dfe17b7f : Traceback (most recent call last):
  File
  "(PATH)\Toolbox\DwgToZippedShapefile.GPServer\extracted\v101\my_toolboxes\UpdateLayer.py",
  line 25, in arcpy.FeatureClassToShapefile_conversion(inFeatures,
  outLocation) File "(PATH)arcgis\server\arcpy\arcpy\conversion.py",
  line 2775, in FeatureClassToShapefile raise e ExecuteError: Failed to
  execute. Parameters are not valid. ERROR 000732: Input Features:
  Dataset Boundary.dwg\Polyline does not exist or is not supported
  Failed to execute (FeatureClassToShapefile). Failed to execute
  (LayerUpdate). Failed to execute (LayerUpdate).

However I don't know why it isn't working as I am passing all the parameters and when running the tool in ArcMap with the exact same parameters it works fine.


texts: Creating and interpolating raster stack from CSV in R

texts:
I have a CSV of biweekly temperature values for the year 2011, at a 1000-m scale in part of the state of Iowa.  The CSV contains ~95000 rows of temperature data, with the first 2 columns as the x and y points and the subsequent 27 columns are the 14-day intervals.  I am trying to convert the CSV into a raster stack at 30-m resolution using interpolation, where each layer of the stack is a biweekly period (27 layers).  The code I am using to convert from a CSV to raster stack also uses a reference raster, which is a raster layer of the correct extent and 30-m resolution.
This code has worked for me before for smaller CSV files.

pts <- read.csv("bimonthly_1000m_temp.csv")

ref = raster("~/reference_raster.tif")
u=crs(ref)

coordinates(pts)=~y+x
proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
pts = spTransform(pts,crs(u))
gridded(pts) = TRUE
r = stack(pts)
projection(r) = crs(u)
r_crop <- crop(r,ref)
extent(r_crop)=extent(ref)
r_crop_resam = resample(r_crop, ref, method='bilinear', filename = "~/interpolated_raster.tif")

But, I get an error at gridded(pts) = TRUE:
suggested tolerance minimum: 0.111429
Error in points2grid(points, tolerance, round) :
  dimension 1 : coordinate intervals are not constant

Any suggestions?  This answer here is in the vein of what I am looking for, though I am trying to interpolate and creating a raster stack.


texts: Summing Attributes based on ID Field using ArcGIS Desktop?

texts:
I'm trying to sum attributes based on if their ID matches the ID of other attributes in the feature class. To explain further, each road has it's own Replacement Cost and Section Number as well as a Super Section ID associated to it. A road can be a part of a Super Section with multiple other roads which will all share the same Super Section ID. What I want to do produce the total replacement cost of every Super Section for each road.

As you can see in the image above, the field Super_Section_Replacement_Cost has the replacement cost for the entire Supersection (the sum of the individual road replacement costs that are apart of the SS).

I was thinking an UpdateCursor might be useful but I'm not sure how I can group the attributes based on a specific field (the Super Section ID).


texts: I want to snap my point feature with Line features. So That I can perform join attribute by location to read attributed from points into lines. Is there any way I can do this task automated. Once points will be snapped/near to lines, I can do rest. Please respond if anyone has an idea.


texts: How to convert from GeoTIFF to PNG/JPG and Web Mercator in Java?

texts:
I'd like to convert raster image from (geo)tiff to an  png/jpg and afterwards if it is necessary ... I'd like to convert the coordinates with help of EPSG-Codes to the OSM/Google Projection(EPSG: 3857/ 900913). All this i'm doing with java programming.
Anyone a good hint... like an java libary??

Thanks a lot.


texts:
My Tilemill exports (MBTiles) tend to fail quite often lately... (Resume Tilemill MBTiles export after crash a bounty is waiting for anybody who helps me there btw ;-) ). We kind of solved the problem (at least a bit) by adding more RAM to our machine...

But still... after an export of a quite small sample area I get the following output:

[3h 57m 23s] Part(1/1) 100.0000% 411.9k/411.9k @ 79/s | 0s left | ✓ 411.6k ■ 39 □ 0 fail 225

So I have 255 failed Tiles. I also have my "log" file cause I started the export adding a --job parameter. In the -failed log file I can see the coordinates of the failed tiles...

Is there a way I can tell Tilemill to redo only those failed tiles?


texts:
I got code to densify my shapefile but my output got empty columns. How to add to ogr2ogr command request to save input columns to output file in following code?

target.write('ogr2ogr.exe -f "esri shapefile" -dialect sqlite -sql "select DissolvePoints(ST_Segmentize(%s.geometry, 1)) as geometry from %s" "%s\\sample.shp" "%s" -explodecollections \n' % (name, name, setupPath, shpin))

I tried to add to sql expression SELECT * , but it doesn't work.


texts: Converting NetCDF WRF weather forecasts from a custom Lambert Conformal Conical projection to WSG84

texts:
I'm having some trouble reading a WRF weather forecast prediction file where the coordinates are in a custom Lambert Conformal Conical projection.
My final goal is to get the weather forecast a a certain location given in a WGS84 projection (lat/lon), which I know for sure is within the NetCDF file.
I'm using a Python script and the NetCDF4 and PyProj libraries.
In the NetCDF file, downloaded from the local weather forecast service, there is the following information about the used projection:
<class 'netCDF4._netCDF4.Variable'>
int32 Lambert_Conformal()
    grid_mapping_name: lambert_conformal_conic
    longitude_of_central_meridian: -14.100000381469727
    standard_parallel: [43. 43.]
    latitude_of_projection_origin: 40.71900177001953
    false_easting: -241.44951
    false_northing: -4.10447
unlimited dimensions:
current shape = ()
filling on, default _FillValue of -2147483647 used

So I try to recreate this using the pyproj library:
lat_0 = variables['Lambert_Conformal'].latitude_of_projection_origin
lon_0 = variables['Lambert_Conformal'].longitude_of_central_meridian
lat_1 = variables['Lambert_Conformal'].standard_parallel[0]
lat_2 = variables['Lambert_Conformal'].standard_parallel[1]
x_0 = variables['Lambert_Conformal'].false_easting
y_0 = variables['Lambert_Conformal'].false_northing

I create a LCC origin projection with PyProj, as well as the destination WSG84 projection:
lcc = Proj(proj='lcc' , lat_1=lat_1, lat_2=lat_2,
           lat_0=lat_0, lon_0=lon_0, x_0=x_0, y_0=y_0,
           a=6370000, b=6370000)
# Resulting projection string: '+proj=lcc +lat_0=40.7190017700195 +lon_0=-14.1000003814697 +lat_1=43 +lat_2=43 +x_0=-241.44951 +y_0=-4.10447 +R=6370000 +units=m +no_defs'

wgs84 = Proj(proj='latlong', datum='WGS84')

Now that I have everything set-up, I perform the actual conversion between projections:
transformer = Transformer.from_proj(proj_from=wgs84, proj_to=lcc)

# Geographical point which I know is contained within the weather forecast.
lat_geo = 42.666981
lon_geo = -8.006540

coords=transformer.transform(lon_geo,lat_geo)

As result, I get:
(497466.92765635485, 234691.1045065646)

which makes no sense to me, as I was expected something contained within the range of the variables in NetCDF file. For reference, the latitude variable in the NetCDF file is between 40.719 and 40.44141.
So, what am I missing here?

Are the lat_0, lat_1, lat_2... parameters correctly loaded? Did I misplaced some of them?
Is there some missing projection information in NetCDF file?
I see that there is actually a single standard parallel (the same, but repeated twice). Does that make sense?

TLDR: How do I go from this custom LCC projection to standard WGS84 coordinates?
More information:

The weather forecast agency says that the global lon/lat coordinates are supossed to be:

GeospatialCoverage:

    Longitude: -11.261331695726343 to -5.101928395629246 degrees_east
    Latitude: 40.421302033295966 to 45.230011083267385 degrees_north

This is a map using their available web tools:


texts: Finding angle to north on image with three control points

texts:
So I am working on some software that geolocates spots on a map taken from a drone. I am having some trouble with the actual geolocation when the map that the program is given is not aligned north.

So what I am trying to do is take the three control points and using them determine at what angle (from the vertical of the image) north is at and then rotate the image by that angle.

What I am currently doing is using GeographicLib (c++ library) to determine the azimuth between two of the control points then the angle in respect to the top of the image of those two and rotating the image by the difference, this has yet to prove correct although i get some numbers close to what they should be if I us the forward azimuth at control Point #2.

I know what the map should look like when aligned as I loaded it into QGIS and rotated it. One other issue is that QGIS seems to have warped it a bit when we overlayed it on a georeferenced image so I do not know if I have to take that into account either.


texts:
I've got a serious problem with the following Grass error:

r.in.ascii -si input=srsaviLC815173.asc output=srsaviLC815173 --o --q
WARNING: closecell: can't move
     C:\Users\bla\Documents\grassdb/utm32/PERMANENT/.tmp/4732.1
     to cell file
     C:\Users\bla\Documents\grassdb/utm32/PERMANENT/fcell/srsaviLC815173

The weird thing is that the same action has been performed already several times with no error. Has anyone got a hint where this comes from?


texts: How to add a GeoJSON to a map with Leaflet

texts:
I need to add a GeoJSON of a shape to a map using Leaflet.
My code:
var map = L.map('map').setView([43.42851395329958, -3.668450798388781], 23);

var geojsonFeature = {
   "type":"Feature",
   "geometry":{"type":"MultiLineString","coordinates":[[[445886.845,4808681.6768],[445887.5215,4808682.5332]]]},
   "properties":{"id":"nudo_153-acometida","x1":445886.8449913647,"y1":4808681.676796015,"x2":445887.52153996995,"y2":4808682.533211502,"longitud":1.091405287449764,"agrupadora":"AGR_AC BT RZ25 P","conductor":"RZ 0,6\/1 kV 3x25+54,6 Alm","ct":"946582","elemento_c":"CABLE","feeder":"33882643","imax":100.0,"nudofin":"acometida","nudoorigen":"nudo_153","numcircuit":3.0,"posicion":"final","r":1.2,"tension":400.0,"tipo":"aereo","tipo_actua":"extension","tipo_eleme":"acometida","tipo_tendi":"posado","tipotramo":"bajada","x":0.1},
   "id":"red_extension.3"
 };
 var myStyle = {
 "color": "yellow",
 "weight": 6,
 "opacity": 1
 };

 L.geoJSON(geojsonFeature, {
   style: myStyle
 }).addTo(map);

 var circle = L.circle([43.42851395329958, -3.668450798388781], {
 color: 'red',
 fillColor: '#f03',
 fillOpacity: 0.5,
 radius: 100
}).addTo(map);

 var wmsLayer = L.tileLayer.wms("http://ovc.catastro.meh.es/Cartografia/WMS/ServidorWMS.aspx?", {
    layers: "Catastro",//nombre de la capa (ver get capabilities)
    format: 'image/jpeg',
  transparent: true,
  maxZoom: 19,
    version: '1.1.1',//wms version (ver get capabilities)
    attribution: "DIRECCION GENERAL DEL CATASTRO"
 });

wmsLayer.addTo(map);

The result is as follows:

Why is the GeoJSON not showing on my map?


texts: Documenting custom Python functions in QGIS using styles?

texts:
This question shows how to document a custom function. However, I'd like to document my function in the same way the rest of the functions are.

So far, I just can manage to create the highlited part of the documentation but I can't find anywhere how to implement the same format used for the rest of the functions.

Is it possible to accomplish this in a custom function?


texts: QGIS2web text labels not visible on the map

texts:
I want to have the text labels visible on my map.
Unfortunately, I don't know where is the error.
In both cases, the labels are based on one category. In one object I can see the label, in other not. When launch the map preview situation is opposite.
Doesn anyone know how to display the text labels in QGIs2web properly?


texts: What are some 3D visualisation tools that successfully work with GIS data?

texts:
There are a large number of 3d modelling tools available (both open and closed source).  I would like to be able to:

Import GIS data into a 3d environment and retain georeferenceing.
Retain the quality and resolution of the source data.
View the data in real time with game like speed.

I have tried a number of tools but they never quite seem to be able to do these three points successfully.  Do any tools exist that can reliably satisfy these requirements?


texts:
I want to save shapefile (polygon) with color in gmt format and in geoghrapic coordinate not UTM. Mmy shapefile is in UTM. I run below as an example how to show in lat lon and with color.

ogr2ogr -f "GMT" Hudson_bounds.gmt  Hudson_bounds.shp


texts: How to make pen style show up correct in QGIS map composer legend?

texts:
As far as I can tell, yes I've updated the legend. See new screenshot. It's like the field width is cutting them, they all begin with one long and then one short line..

I have a several different line shapes files that represent buss lines. When I set "Pen style" to different styles I see the changes on the map but it doesn't appear correct in the legend in QGIS map composer.

Fix or work around needed! Using QGIS 2.0 on Mac.


texts: Projection issue - layers misaligned despite same coordinate system (British National Grid) in ArcMap

texts: I'm having projection issues with a shape file containing habitat information not lining up with my other layers. This is particularly evident when comparing the shape file to the aerial imagery layer I'm using, which is downloaded from Digimap (see image of water bodies). All layers are using the same coordinate system, British National Grid. I have tried setting the coordinate system then importing the files in different orders with the same result. I'm working in ArcMap 10.7 but I've also opened the files in QGIS with the same issue. I'm a GIS beginner, is there anything else I should check or do to resolve this issue?


texts: geopandas geometry.within problem ValueError: Geometry column cannot contain mutiple geometry types

texts:
I'm trying to do an overlap of some lines with polygons with geopandas.

For a reason I don't understand, gpd.overlay(grid, hel, how='intersection') is not made for lines and polygons, so I'm trying to do it like this:

final = grid[grid.geometry.within(hel)]
  final.to_file('output.shp')

but I get this error: ValueError: Geometry column cannot contain mutiple geometry types when writing to file. I don't see how this processing can give any different geometry type rather than line.

If I the overlay in QGIS it works perfectly.

Any ideas?

EDIT
As pointed in the comments, this error comes from to_file command, caused by trying to write an empty GDF. This error can be reproduced by :

g = gpd.GeoDataFrame()
g['geometry'] = None
g.to_file('testburro.shp')

I will open an issue about this, the error message is totally misleading.

Back to my problem, this code produces the desired output:

with fiona.open(filename1, 'r') as layer1:
with fiona.open(file_name, 'r') as layer2:
    index = rtree.index.Index()
    for feat1 in layer1:
        fid = int(feat1['id'])
        if feat1['geometry']:
            geom1 = shape(feat1['geometry'])
            index.insert(fid, geom1.bounds)
    drv = layer1.driver
    crs = layer1.crs
    schema = {"geometry": "LineString", "properties": {"id": "int"}}
    output_name = 'output.shp'
    with fiona.open(output_name, 'w', driver=drv, schema=schema, crs=crs) as dst_layer:
        for feat2 in layer2:
            geom2 = shape(feat2['geometry'])
            for fid in list(index.intersection(geom2.bounds)):
                feat1 = layer1[fid]
                geom1 = shape(feat1['geometry'])
                if geom1.intersects(geom2):
                    geom3 = geom1.intersection(geom2)
                    dst_layer.write({'geometry': mapping(geom3), "properties": {"id": feat1}})
                    print('{} intersects {}'.format(feat2['id'], feat1['id']))

So I don't think is a problem of the dataset.


texts: How to delete rings automatically using QGIS/Grass

texts:
I'm cleaning up a (very) messy shapefile in QGIS/Grass.

I now have several large polygons with numerous (100s) of small rings.

Is there a quick and easy way to remove/ fill all these small rings. Selecting and deleting each one manually is not feasible.


texts: How to interact with other layers in a QGIS Python init function?

texts:
I am using QGis and making a Python script that sets the Python Init function to a custom UI file. I have a function with the following signature:

def formOpen(dialog, layer, feature):
    ...

I can access the given layer and it's contents - this works fine. But what should I do in order to access other layers?

My use case is the following: Whenever a given feature F of layer A is modified, I want to calculate the intersection of feature F with all geometry objects of layer B and show it to user. This information is later saved to the database.

I understand that in order to do that, I should be able to access the interface object. In my plugin I save the reference to the interface by:

self.iface = iface.

However, the Python init function is out of my plugin. Should I put it inside my plugin? Or use global variables?


texts: Finding missing dependencies in distributed ArcGIS Python toolbox?

texts:
I am developing a Python Toolbox that uses a few external libraries (ogr/gdal, matplotlib, shapely, reportlab among others) which files I have directly copied in the same folder as my toolbox .py, as recommended in Esri's A simple approach for including 3rd party Python libraries with your scripts. It works just fine on my computer, but when I run this tool on another computer, I get an ImportError. It's very likely that the libraries are the reason for that error.

<type 'exceptions.ImportError'>: DLL load failed: The specified procedure could not be found.

Moreover, I don't get the ImportError the same way from ArcGIS10.0 or ArcGIS10.1 (10.0: only get the error without any additional message; 10.1: get the error where it worked fine in the previous version), which worries me a bit because this really looks like a singular configuration dependent problem.

How do I find the missing dependencies in my current distribution and ensure that the tool which I distribute will work on every computer with a fresh ArcGIS installation?

I have tried what I consider the usual ways to solve that problem:

Run the script outside ArcMap to debug it but I did not get much more information.
Use Dependency Walker to find the missing DLLs but when I try to profile ArcMap, Dependency Walker stops running.


texts:
I am trying to mosaic some raster-datasets (~15) using rasterio.merge.merge in python.
I defined a method to calculate the mean but have trouble to take the no-data-values into account.
def nanmean(old_data, new_data, old_nodata = -10000, new_nodata = -10000, index=None, roff=None, coff=None):
    old = np.where(old_data == old_nodata, np.nan, old_data)
    new = np.where(new_data == new_nodata, np.nan, new_data)
    old_data[:] = np.nanmean([old, new], axis = 0)

mosaic, out_trans = merge(src_files_to_mosaic, method= nanmean)

That is my function, if I use first/last or min/max as method everything seems to work perfectly fine.


texts: Making several SVG icons using one layer in QGIS?

texts:
I need to make several different SVG icons in the same size.
I am using a list with numbered SVG from 1 to 262:

The obvious way is to make a different layer for each SVG.

I have tried to change the SVG icon from the Layer Styling bar, but whatever I did I only get the same SVG over and over again

Is there a way to use different SVG with one layer?


texts:
We are trying to use proj4 to project a local coordinate system to EPSG4326. We believe we have most of the string created correctly in QGIS and as follows;

+proj=ortho +datum=WGS84 +lat_0=40.52741251 +lon_0=-117.12455304 +x_0=372979.646 +y_0=103881.838 +units=us-ft no_defs

We are missing azimuth and scale. I believe we can use +k or _k_0 for scale but cannot determine how azimuth is defined. Are these incorporated into such a projection in proj4?
We have also run the ESRI equivalent prj through GDAL srsinfo with the following output:

PROJ.4 : ERROR 6: No translation for LOCAL to PROJ.4 format is known.

Local system params in ESRI.prj;

Projection: Local
Scale_Factor:   1.00024937
Azimuth:    -0.31500000


texts: Geoserver - Layer creation via data import

texts:
We have about 400 sites (materialized by coordinates in the spreadsheet file to be imported) and several fields of information attached to these sites (coordinates, site ID, name, comments, category, etc). This data have been collected during a site survey and now need to be presented on a map, with details shown when a user clicks on the site spot appearing on the map.

Is there any particular workflow that one would recommend to successfully complete this task with Geoserver and avoid a manual entry of all the data already available on a well structured Excel spreadsheet?


texts: Controlling for 1st order effect when testing for spatial association between points of two different patterns

texts:
Background

I have a follow-up question relative to a question I asked long time ago link to earlier question. I am asking a new question since this is something new, even though it builds on the earlier thread.

Aim

In my earlier thread, I was interested in formally (statistically) assessing if the points of a pattern A tend to be close to the point of another pattern (B). I got an extremely useful answer that was suggesting to perform a permutation-based procedure (notably in R, which sounds perfect to me) to set the observed average minimum distance against the distribution of average minimum distances generated under the Null Hypth of no "association" between the two patterns.

Now, I wanted to know if there is a way for controlling for a 1-st order effect, that is for the influence of an underlying spatial covariate. In other words, points of pattern A and pattern B could be "associated" because they tend to "sit" on particular patches of land (e.g., type of soil, particular elevation values, slopes, etc).

Question

How can I control for the effect of a spatial covariate?

Is it statistically sound to keep (say) pattern B fixed and generate (say) 199 random A points, each time calculating the average minimum distance between the points of the 2 patterns, and then work out a p-value? The generation of random A points would be done using the spatial covariate as probabilty density function.


texts: How to convert Landsat L1G raster to Geotiff?

texts: Is it possible to convert a Landsat L1G raster file like L71001070_07020000728_B30.L1G to Geotiff using free software?


texts: Getting text name of current project CRS in QGIS Layout

texts:
I'm want to display the current project projection in a QGIS layout.

Right now I can get QGIS to display the current EPSG number like below:

But what I want to do is get QGIS layout to print the ESPG name as well like below with an expression:


texts: how to load the QGIS openlayers_plugin in standalone python script (not the QGIS embedded python console)

texts:
if I try to load the openlayers plugin in a standalone application
it always fails with False, although paths are set, and plugin shows available.

(tested on Ubuntu 14.04, QGIS 2.4, python-2.7)

#!/usr/bin/python
# -*- coding: utf-8 -*-

import sys, os

import qgis
import qgis.gui, qgis.utils
from qgis.core import *
from PyQt4 import QtGui

# open an app
app = QtGui.QApplication(sys.argv)
# supply path to where is your qgis installed
QgsApplication.setPrefixPath(u'/usr', True)
# load providers
QgsApplication.initQgis()

# set up the plugins
qgis.utils.plugin_paths = \
    ['/usr/share/qgis/python/plugins',
    os.path.expanduser('~/.qgis2/python/plugins'),]
qgis.utils.updateAvailablePlugins()
print qgis.utils.available_plugins
print "...load:", qgis.utils.loadPlugin(u'openlayers_plugin')
print "...start:", qgis.utils.startPlugin(u'openlayers_plugin')
print "active:\n",qgis.utils.active_plugins

canvas = qgis.gui.QgsMapCanvas()
canvas.show()

what is the issue here ?
post How to fetch openlayers layers from pyqgis?
refers to the embedded python console, where the plugin was magically available before.

thanks!


texts:
I want to lable  feature class with related table.One of the solution for this goal is writing script in the Expression of label tab ( How To: Label a related table). But i don't want to scripting in Expression.Then it's not as simple as the Esri solution.

My code is working properly but i could not iterate "lblclass.expression" for the row1 and  the related table ( row1) is not iterate and just show one value for all labels.How to iterate related table value for labeling ?

import arcpy
table = "table"
listtlabel = []
mxd = arcpy.mapping.MapDocument("CURRENT")
layer = arcpy.mapping.ListLayers(mxd, "Point")[0]
fielddelimi1 = arcpy.AddFieldDelimiters(table, "Code")
sql_exp1 = "{0} = {1}".format(fielddelimi1,309)

if layer.supports("LABELCLASSES"):
    for lblclass in layer.labelClasses:
        with arcpy.da.SearchCursor(table,["Code1","Row"],sql_exp1) as cur1:
                for row1 in cur1:
                    listtlabel.append(row1[0])

                lblclass.expression = '[Code]' + '&'+ str(row1[1])
                str1=','.join(str(e) for e in listtlabel)
                print str1
                with arcpy.da.SearchCursor(layer,"Code2") as cur:
                    for row in cur:
                        if row[0] in listtlabel:
                             lblclass.SQLQuery = '"Code2"' + "IN" + "({0})".format(str1)
lblclass.showClassLabels = True
layer.showLabels = True
arcpy.RefreshActiveView()
del mxd


texts: How to reclassify a raster based on quantile and without no data in QGIS?

texts:
I want to reclassify quantile a raster layer without negative values (nodata values) so I used Reclassify a raster file with quantiles and How to filter no data value with GDAL? in below code:

import numpy as np
from osgeo import gdal, gdal_array

# open the dataset and retrieve raster data as an array
dataset = gdal.Open("path/to/file.tif")
array = dataset.ReadAsArray()
lst = []
for v in array:
    for element in v:
        if element > 0:
            lst.append(element)

# create an array of zeros the same shape as the input array
output = np.zeros_like(lst).astype(np.uint8)

# use the numpy percentile function to calculate percentile thresholds
percentile_80 = np.percentile(lst, 80)
print percentile_80
percentile_60 = np.percentile(lst, 60)
print percentile_60
percentile_40 = np.percentile(lst, 40)
print percentile_40
percentile_20 = np.percentile(lst, 20)
print percentile_20
percentile_0 = np.percentile(lst, 0)
print percentile_0

output = np.where((lst > percentile_0), 1, output)
output = np.where((lst > percentile_20), 2, output)
output = np.where((lst > percentile_40), 3, output)
output = np.where((lst > percentile_60), 4, output)
output = np.where((lst > percentile_80), 5, output)

outname = "/path/to/newfile.tif"
gdal_array.SaveArray(output, outname, "gtiff", prototype=dataset)

But I received below error:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/home/nikan/Untitled-2.py", line 35, in <module>
    gdal_array.SaveArray(output, outname, "gtiff", prototype=dataset)
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal_array.py", line 239, in SaveArray
    return driver.CreateCopy( filename, OpenArray(src_array,prototype) )
  File "/usr/lib/python2.7/dist-packages/osgeo/gdal.py", line 1538, in CreateCopy
    return _gdal.Driver_CreateCopy(self, *args, **kwargs)
ValueError: Received a NULL pointer.


texts: How to use ST_ClusterWithin with geography data type and measurement units?

texts:
I want to clustering locations use ST_ClusterWithin function:

geometry[] ST_ClusterWithin(geometry set g, float8 distance);

My table has column location of geography(Point,4326) type and I don't understand how to convert it correctly to geometry type with ability to specify distance in meters.

The following query can convert geography to geometry but allow to specify distance in degrees that is not acceptable for me:

select ST_AsText(unnest(ST_ClusterWithin(location::geometry, 0.0001)))
from places

So, how to use this function with geography data types?


texts: While running SAGA 7.30 version in QGIS 3.16, they show an error about that version is not supporting and checking about a new version for up-gradation that might solve the issue but 7.30 version is the newest. How can I fix this problem?


texts: How to label only some places instead of all?

texts:
How do I turn of some place names, or rather turn on some places names but not the rest.  How do I turn of or hide that?

I also get a red cross under the place names which I don't like to see.
How do I turn of or hide that?

The place names is in a line layer (using labelling), and there are also many numbers I  don't need to see.  Some of them are Height above mean sea level (AMSL), but not all.

QGIS Lisboa versjon 1.8.0


texts: Exporting selected features in dxf using PyQGIS

texts:
This code should create a new DXF file for every new features contained in an attribute of a shapefile.
layer = QgsProject.instance().mapLayersByName('LAYER NAME')[0]
filename = 'FILE NAME'

objects = []

for l in layer.getFeatures():
    objects.append(l["ATTRIBUTE NAME"])

object = set(objects)

print(objects)
print(object)

for o in objects:
    layer.selectByExpression("\"type\"= {}".format(o))
    QgsVectorFileWriter.writeAsVectorFormat(layer, filename + o + '.dxf', 'utf-8', layer.crs(), driverName='DXF', onlySelected=True, skipAttributeCreation=True)

But the following error occurs (the dxf file is created on the computer, but it doesn’t contain any information) :
(5, "Creation of field ID failed (OGR error: DXF layer does not support arbitrary field creation, field 'ID' not created.)")
As mentioned in this post Export DXF from QGIS layer using PyQGIS?, if skipAttributeCreation=True is set, no error occurs, the dxf files are created but there’s no information…
(0, '')


texts: QGIS GRASS Processing Not Working on MacOS?

texts:
Recently upgraded to QGIS3 (3.0.2) for MacOS (High Sierra) and encountered GRASS (7.2.2) issues. Tried adding a custom variable and now when running GRASS algorithms (e.g. r.reclass) I get the following error:

WARNING: Default locale settings are missing. GRASS running with C locale.
Starting GRASS GIS...
ERROR: Path '//' doesn't exist
Default locale not found, using UTF-8
Exiting...

The GRASS7 Path in Settings --> Options --> Processing seems to be correct. What am I missing?


texts: Attribute table not updating after QGIS field calculator used?

texts:
Using QGIS 2.8.1 i am trying to update an attribute table using the field calculator.  its a vector layer of vegetation communities.  i am trying to update the area fields in the table.  The output example in the calculator shows a correct calculation, but the areas do not update and stay as 0.0000.  The field is "real" 10 wide with a precision of 4.

Any ideas on why the table does not update?


texts: How to add spatial reference ESRI:102100 as a custom CRS to Geoserver?

texts:
I've been struggling with this for over a week. I'm using the ArcGIS Javascript API to create a basic map application, and I'm trying to add a WMS layer from my geoserver to this application. The problem is that the basic map app I'm working with is using an ESRI basemap with spatial reference ESRI:102100. My wms layer is using spatial reference EPSG:4326. See my other question for more details about this: ArcGIS Javascript API - How to define spatial reference for WMS layer?

When I run my application I get this error:

<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.3.0" xsi:schemaLocation="http://www.opengis.net/ogc http://localhost:8040/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">
<ServiceException code="InvalidCRS">
Error occurred decoding the espg code urn:x-ogc:def:crs:EPSG:102100 No code "EPSG:102100" from authority "European Petroleum Survey Group" found for object of type "IdentifiedObject".
</ServiceException>
</ServiceExceptionReport>

I tried to add ESRI:102100 as a custom CRS to my Geoserver (by editing the epsg.properties file), but it adds it as EPSG:102100 which is not a valid epsg code according to the error message I'm getting. Anybody know how to get Geoserver to recognize ESRI:102100? Or, anybody know how to properly define a wms layer from a geoserver in ArcGIS Javascript API?


texts:
I am trying to import some line data representing utilities into GE Pro. However they get a really high offset (pic 1). I've tried changing the height in the "style template" menu that you can choose when you import the shp data. "Hoogte" is height in my language and I've put it to 0. However as you can see in pic 1 they keep having a really high offset. I've tried setting "Hoogte" to all sorts of different values to no avail. In pic 3 I have another menu from you can navigate to by right clicking the layer in GE pro. Here I am able to successfully change the height, but it messes up the colour scheme and I can't put it back to the required colours.
TLDR: I want to put this shp data at ground level and not have a 50m offset.


texts:
I need help with ArcGIS API for JavaScript, I need to use private modules in code for my tools, but I can't figure out how use them in my app like on click launch this module.

I'm lost in DoJo...

The map.js:

     /**
 * This file contain :
   * - Map init
   * - Principal layers to add
   * - Reference to the tools files
   * - Style of the principal page (map, navbar...
  **/

require({
    packages: [{ name: "src", location: "/src", type: "text/javascript" }]
},[
"esri/map",
"esri/toolbars/draw",
"esri/graphic",
"dojo/parser",
"esri/symbols/SimpleFillSymbol",

"dijit/registry",
"esri/geometry/Geometry",
"dijit/layout/BorderContainer",
"dijit/layout/ContentPane",
"dijit/form/Button",
"dijit/WidgetSet",
"src/modules/zoningSelection",
"dojo/domReady!"],
function(Map, Draw, Graphic, parser, registry, dijit, zoningSelection){
   parser.parse();

   map = new Map(map,{
      basemap:"hybrid",
      center:[-1.486075,46.199997],
      zoom:15
   });
   console.log(map);
    var selectionTool = dijit.byId("select");
    selectionTool.on("click",startSelect);

    var validationTool = dijit.byId("end");
    validationTool.on("click", endSelect);

    var deleteTool = dijit.byId("del");
    deleteTool.on("click", delSelect);

   function startSelect(){

        var select = new zoningSelection(map, Draw.POLYGON);
        console.log(select);
    };

    function endSelect (){
        alert("Validé");
    }
    function delSelect(){
        alert("Suppr");
    }
});

My Module to load in the first code:

require(["esri/map",
    "esri/toolbars/draw",
    "esri/graphic",

    "esri/symbols/SimpleFillSymbol",

    "dojo/parser", "dijit/registry",
    "esri/layers/GraphicsLayer",
    "esri/graphicsUtils",
    "esri/geometry/Polygon",
    "esri/geometry/Geometry",
    "esri/dijit/Popup",
    "dijit/layout/BorderContainer", "dijit/layout/ContentPane",
    "dijit/form/Button", "dijit/WidgetSet",
    //"dojo/_base/array",
    "dojo/domReady!"],function (map, drawType, Draw, GraphicsLayer, Graphic, SimpleFillSymbol, graphicsUtils, dijit) {
        console.log("séléction lancée ou non ");
        console.log(map);
        console.log(drawType);
    var _drawTool = new Draw(map);
    var _drawingLayer = new GraphicsLayer();
    console.log(_drawingLayer);

    map.addLayer(_drawingLayer);

    _drawTool.setRespectDrawingVertexOrder(true);
    _drawTool.on("draw-complete", _onDrawComplete);

    activate(drawType);
    //map.disableMapNavigation();

    var _onDrawComplete = function(evt) {
        var symbol = new SimpleFillSymbol();
        map.enableMapNavigation();

        var geom = evt.geometry;
        //create object from geometry and symbol
        var graphic = new Graphic(geom, symbol);
        _drawingLayer.add(graphic);

        var points = graphicsUtils.getGeometries(_drawingLayer.graphics)[0].rings[0];

        require(["dojo/_base/array"], function(array){
            array.forEach(points, function(point, i){
                console.log("point #" + i + " : x=" + point[0] + ", y=" + point[1]);
            });
        });
    };
});

How can I use them together?


texts: How to find points that lay on polygons in which another point layer exists by using FME 2021?

texts: I have 3 point feature classes (A,B and C) and 1 polygon feature class. I am trying to extract only features from feature class A that lay inside polygons which already containg features from feature class B or C. I tried "PointOnAreaOverlayer" and "SpatialRelator" Transformers but for some reason I dont get desired results. Any hint how this can be done in FME 2021 ?


texts: Saving a form .ui file with relative path?

texts:
My company uses QGIS 2.18 (soon to upgrade to 3.4) and I've designed a form to make data entry for our survey field staff a little less hateful.

The .ui file is saved on a network drive, and as long as our field staff have access to the network, the form appears when they're doing data entry. However, sometimes they're working pretty remotely, with no access to a network (or the internet at all) and so the form no longer appears.

Is there a way to set up the .ui file so that it can be copied along with the other shapefile files to a users C:\ drive and have a relative path setting so that the .ui file will always be referenced no matter what folder the shapefile/ .ui file is saved in?

The workaround that I've thought about trying is to copy the .ui file to each users (around 15 ppl) computer and set up the shapefile .ui for each individual person, however this seems a bit cumbersome, will be prone to error, and just seems a bit inelegant.


texts: Using python to do a calculation on a set of numbers in a text field separated by commas

texts:
Using ArcGIS ArcMap 10.6:
I have a text field in an attribute table that has comma separated values (example: 12,24, 64, 87,29 OR 13, 45, 32,21,34, 45). For this example, I will name this field VALUES. Please note that VALUES includes data collected from field work so the spacing between commas will be varied. My goal is to be able to run a calculation on each number before the comma, and to get the sum of that number in a new field called Calculated_Values.
Though, I am currently using excel, as I cannot figure out how to do this with python. My current workflow is as followed:

Export the ID and VALUES fields in the attribute table into a .csv file using Table to Table in ArcMap 10.6
Open the .csv in Excel and run the 'text to columns' function under the Data ribbon
Separate the VALUES field into multiple columns based on commas
Run this calculation: =ROUND(SQRT(POWER(A2,2)+POWER(B2,2)+POWER(C2,2)+POWER(D2,2)+POWER(E2,2)+POWER(F2,2)+POWER(G2,2)+POWER(H2,2)),0) in Column I
Copy Column I values and ID columns into a new .csv file
Save and close the new .csv file and join the shapefile to the new .csv file based on ID
Export the newly joined shapefile into a new file and delete the duplicate ID field in the new shapefile

I am wondering if there is a way to automate this Python/ArcPy so I can create a custom script tool to do this without leaving the ArcGIS environment.
Is there a way to separate values in a text field based on commas into different columns, run the above calculation, and have the values from that calculation entered in a new field named VALUES_CALCULATED?
I know I have to use a cursor but don't know what type of for loop to run in it.


texts: Can Hyperlink be based on two fields for path and file name separately?

texts:
I need to create a hyperlink field. I have one field with the file name, and another field that contains only part of the file path (the real path is unknown). The table contains 14,000 records with 14,000 different DWG file names.

Is it possible to create a Hyperlink field in a table or feature class that will start  to search automatically with windows browser for the file at a specific folder according to the field with the estimated Path?  i.e. one field with the file name to look for, and another field with part of the path to look at.


texts: Obtaining trial version of ArcGIS Publisher?

texts:
I would like to save my MXD files as PMFs. I don't have a license of ArcGIS, only the 60-day trial.

Is there a similar trial for the Publisher?

I understand that's the only way to save PMF files (through the Publisher).


texts: Can not install Microsoft SQL Server 2012 Native Client on windows 10

texts:
I am trying to compile QGIS project myself.
I am following the instructions below, which seems to be official guide but outdated.

https://htmlpreview.github.io/?https://github.com/qgis/QGIS/blob/master/doc/INSTALL.html#toc15

When I install spatialite library through osgeo4w-setup-x86 : pyspatialite, spatialite, spatialite4dll.
During the install, Microsoft SQL Server 2012 Native Client is required to install but can't proceed because SQL Server 2012 is not supported on the oparting system (Windows 10/7 )

The libraries seems to be installed all right but SQL Server 2012 Native Client.

Is there any workaround to install SQL Server 2012 Native Client.

Eventually I can't do cmake configure to create output.