texts: Defining a user defined CRS to a layer in QGIS modeler and then reprojecting it

texts:
How do I get the QGIS modeler to recognize and use my USER DEFINED CRS in proj4 format?
I have a user defined CRS with the authid USER:100024.  Its parameters in proj4 are "+proj=utm +zone=18 +ellps=clrk66 +datum=WGS84 +towgs84=-11.0584,163.8703,190.541,0,0,0,0 +units=m +no_defs +wktext".  Manually I can set my user defined CRS on a given layer (right click layer in TOC, 'Set layer CRS' to 100024) and then save to EPSG:2959 (NAD83 UTM18N) and get a valid return of point coordinates.  This works fine.
The problem I'm running into is in trying to replicate this behaviour with a QGIS model.
From the modeler template I first use the "define current projection" on the source layer (where I select USER:100024 from the CRS window).
Second, I run a "reproject layer" to EPSG:2959 on the defined projection from step 1.
But this does not work, the coordinates area about 4m off what they ought to be.  It seems to ignore the defined proj4 layer paramters (unlike when I do it manually).  In fact, when I verify the properties on the source layer it defaults to WGS84 (which sort of makes sense because it's not an official EPSG code) instead of 100024.


texts: Qgis2web filtering clusters with Leaflet output

texts: I'm using the Leaflet output option on qgis2web. However, when I select the 'Cluster' option for my exported layers, the 'Attribute filters' that are created don't appear to update the displayed data / symbols. Have I missed a setting/option or is this not possible with a qgis2web export to Leaflet?


texts: QGIS 2.10.1 (64bit) rotation field setting missing

texts: I cannot find the the rotation field setting under the properties-layer settings-advanced (for styles) menu in the new release of QGIS (2.10.1 64bit). The only setting what is available there is the "symbol levels". Any idea where can I find this function?


texts: WGS72 to NAD83- how to find correct transformation

texts: I am trying to get this raster file to convert from WGS72 to NAD83, but ArcMap10 doesn't appear to have the correct transformation listed.  Where can I find this transformation, and/or how do I import ones that are not already listed as available?


texts: Unique name segment creation in ModelBuilder?

texts:
Is it possible for us to give each output of a process a unique name segment (inline name). So for example I have a process that dissolves, then eliminates and buffers a single input file. I need to save the itermediate steps into a gdb.

So we have [TownNameX_slope]--Dissolve[\townDissolves.gdb\TownNameX_dissolve]--eliminate[townEliminates.gdb\TownNameX_eliminate] etc. then TownNameY etc.

So in a batch process I have to keep entering TownNameX for each output type. I know in iterators you can use %n but this is for just the output of one tool. Can we extend this to all the output where

%n=townnamex and then for dissolve the output is %n_Dissolve etc

I guess the other option is to run it in python in a "for" loop but I would prefer to stay in model builder as the reporting is easier.

This is similar to the question in Using ArcGIS ModelBuilder to perform In-line variable substitution for input data path?


texts: Recovering Latitude and Longitude from Albers to Compute Distance

texts:
I have a shapefile of U.S. counties in USA Contiguous Albers Equal Area Conic.  I want to do distance calculations.  For example, distances from one county to another.

I do not know how to interpret the coordinates.  In another thread, I understand that the distance is in meters.  Is that correct?

The ultimate goal is to use the distance as a part of a statistical procedure.  Specifically, I am using it in a regression discontinuity design, in which I restrict the sample to counties that are a given distance from a particular line in the U.S. (e.g., the Central Time Zone).  So, I'd like the distance to be fairly accurate because I am restricting the sample to be say 20 miles on either side of this line, and the distance should be accurate from a symmetrical perspective.

How may I recover a latitude and longitude?  I wanted to recover them so that I can apply the distance formula for coordinates by hand to get the distances in miles in a way that I can interpret.  I've come across projections and UTM when searching, but I do not know what that means.

I use ArcMap in ArcGIS 10.1.  Thank you very much for any help!


texts: What causes v.net connect "Failed to find intersection segment" and how to fix it?

texts:
In GRASS 7, when running v.net connect, I get the ERROR: Failed to find intersection segment and no output is generated. The only mention of this error Google found is in the source of connect.c

What can cause this error and how can I fix my input data so that I can build a network?


texts: How to Export a Personal Geodatabase to ArcSDE

texts: I want to copy/export my personal geodatabase which consists of feature datasets, relationship class etc to an Enterprise geodatabase. What would be the best approach?


texts:
I am trying to extract the html code out of an ArcGIS Online map from Python.
from arcgis.gis import GIS
import pandas as pd

g = GIS("https://www.arcgis.com",
"name", "password")

p = {
 'type':"CSV",
 'title': '2020#12',
 'snippet':'hwa',
 'description':'test ' + \
 'dataframe object to a GIS item',
 'tags': 'python',

 }

i=g.content.add (item_properties=p,

data="C:/Users/Name/PycharmProjects/colors/Map/New.csv"
             )
q=i.publish ()
print (q)
search =
g.content.search(query="owner:name",
                               item_type="Feature Layer")
print(search)

map=g.map (location="United States")
item=search [0]
print (item)
map.add_layer(item)

props={'title': '2020#156',
 'snippet':'hwa',
 'tags': 'python',
 "renderer":"autocast",
 "field_name": "Threat"}

map.save (props)

map.embed ()

I get an error for called embed ().
AttributeError: 'MapView' object has no attribute
'_preview_image_display_handler'

I have tried the export_to_html () method. It throws no error but the html code is invalid when I try and embed it in google sites.


texts: Error with Gdal 1.10 install on Ubuntu 13.10

texts:
I've installed GDAL 1.9 then upgrade to 1.10 on Ubuntu 13.10 and then now it's error (gdalwarp can work but gdal_translate is fail)

gdal_translate
ERROR 1: gdal_translate was compiled against GDAL 1.9 but current library version is 1.10

How can I fix this problem?

Thanks


texts:
I'm have added a map from a WMS-server in Google Maps v2 for Android using the code provided from this site: http://www.azavea.com/blogs/labs/2013/01/wms-on-android/ . The tiles are downloaded and displayed but something is not working properly. The map looks blurry and the detail level is low.

I have tried fiddling around with the bbox parameters and if I increase the value for each point multiplying by 4, I can see a greater resolution and more details in the map but the tiles are way off.

Do I have to add some parameters to the wmsurl like resolution or something?


texts: Label only the largest polygon of a multi-part polygon layer in QGIS

texts: I am exploring the new Label engine of QGIS 2.12.2 and I found it a great improvement. However, I discover it is not yet up to ArcGIS "standard" when dealing with multipart polygons. As a workaround (unless I am missing some hidden property), is there a way, through an expression to be passed to the label engine, to label only the largest feature (in terms of area or perimeter) of a multipart polygon?


texts: Extracting point data / WMS out of (browser) geoviewer

texts:
How do I extract point data out of a browser-based map?

In particular, I need the park and ride points from http://wegenenverkeer.be/carpoolparkings, but I cannot find their WMS or shapefiles, anywhere. I can load the 'element' and see a script which loads the points, but don't know how to proceed with this code.

It can be WMS, WFS, anything I can load into my ArcMap file.

So the page gives the points in this script, can I use this in arcmap in some way?

{"projection":"EPSG:4326","attributes":{"name":"\u003Ch4\u003E\u003Ca href=\u0022\/carpoolparkings\/oostkamp\u0022\u003EOostkamp\u003C\/a\u003E\u003C\/h4\u003E","description":"Afrit E40"},"wkt":"POINT (3.2319545170254 51.146994290477)"},{"projection":"EPSG:4326","attributes":{"name":"\u003Ch4\u003E\u003Ca href=\u0022\/carpoolparkings\/oostkamp-2\u0022\u003EOostkamp 2\u003C\/a\u003E\u003C\/h4\u003E","description":"N50g \/ A10 (E40) - uitrit 9"},"wkt":"POINT (3.2474898713709 51.14465742022)"},


texts:
I'm trying to create an action in QGIS that open this link

  htt p://www.gulliver.it/itinerario-foto/14451/#!/uploaded/2014/large/20140929123749.jpg||grid2

but the character | is converted into %7C preventing the page beeing open correctly, in Chrome, indeed the result of the action [% "web-link field" %] is:

  htt p://www.gulliver.it/itinerario-foto/14451/#!/uploaded/2014/large/20140929123749.jpg%7C%7Cgrid2

I defined action as you can see in the picture

Here there are the webpages as they appear in my browser:

Rewriting the link in the field and the action like:  [% "truncated web-link field" %]||grid2 is the same.

How is it possible to make the page open correctly?

Alternatively I could open the webpage internally but, following this instructions

  ht tps://docs.qgis.org/2.2/en/docs/training_manual/create_vector_data/actions.html#hard-fa-open-a-webpage-directly-in-qgis
  ,

I haven't any results (nothing happens clicking on elements).

Any solutions?


texts: GDALOpenEx() doesn't works anymore with iOS 12

texts:
I'm using GDAL in a iOS application. Everything was fine since the iOS 12 SDK update. Now, when I try to open a database (.sqlite file here) with this line

GDALOpenEx(dataSourceConnectionString.c_str(), openFlags, NULL, NULL, NULL));

dataSourceConnectionString is equivalent to

/var/mobile/Containers/Data/Application/E2E5852A-F597-40F8-86FC-7BAFA298CCD2/Documents/data/database.sqlite

The path is correct and the file exists. The error number raised by GDALOpenEx is 4, which means the file is not recognized as a supported file format.

It works perfectly well if the app is running on iPad under iOS 11, but it doesn't anymore if it's iOS 12.

Has anyone experienced similar issues with GDAL/sqlite3 and iOS 12 update ?


texts:
I am having some problem with Gdal.open with my list of raster. I use Python 3 with GDAL 2.2.

This is my script:

import os
import glob
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage
import pandas as pd
import datetime as dt

inDir = 'C:/Doy2000/'
os.chdir(inDir)
#repertoir sorti

outDir = os.path.normpath(os.path.split(inDir)[0] + os.sep + 'output') + '\\'
if not os.path.exists(outDir): os.makedirs(outDir)

# Create list of MODIS rasters for surface reflectance
modis_EVI = glob.glob('MOD13Q1.006__250m_16_days_EVI_**.tif')

evi = gdal.Open(modis_EVI[0])

I don't know what I can do. When I do print(modis_EVI ) and print(evi) I had :

print(modis_EVI)
['MOD13Q1.006__250m_16_days_EVI_doy2000209_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000225_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000241_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000257_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000273_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000289_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000305_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000321_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000337_aid0001.tif', 'MOD13Q1.006__250m_16_days_EVI_doy2000353_aid0001.tif']

print(modis_EVI[0])
MOD13Q1.006__250m_16_days_EVI_doy2000209_aid0001.tif
print(evi)
None

What can I do for this error? How can I control this error?


texts: High resolution, printable alternative to OpenLayers plugin for QGIS?

texts:
Are there any alternatives to OpenLayers plugin to have Google / Bing maps for QGIS projects?

I need an high quality cartography (both aerial and street view) but OpenLayers plugin does not allow printing of these.

My configuration is:

QGIS 2.6.1
Windows 7 64-bit
OpenLayers version 1.3.5


texts: Network Analyst resolved route can't display properly on printed map

texts:
I have looked on the Internet but there are no updated answers. I am using ArcMap 10.2. I am trying to print the route resolved by network analyst but there is no way I can display overlapped generated route when printed.

This is an example where the encircled section is passed through three times. Please help with how I can show them separately.


texts: Cost Path tool in script randomly throws error

texts:
I have a script tool that is creating a least cost path for sets of points using a Cursor. When applied, the script runs great, except every so often the tool is stopped because there is an error in the Cost Path tool. If I close out of the tool and run the last few steps manually on the points everything works fine and I can start the tool on the next row.
Here is the cursor in the script:

rows = arcpy.SearchCursor(setExtract)
for row in rows:
    if row.setX == row.upX and row.setY == row.upY:
        arcpy.AddMessage("Points for trawl " + str(uniqueID) + " have identical coordinates, the trawl cannot be drawn.")
    else:
        avgDepth = row.avg_depth
        uniqueID = row.uniqueID
        arcpy.MakeFeatureLayer_management(setExtract, "setpoint", "uniqueID = " + str(uniqueID))
        arcpy.MakeFeatureLayer_management(upExtract, "uppoint", "uniqueID = " + str(uniqueID))
        outRas = Abs(Raster(clipBathy) - avgDepth)
        outRas.save(tempRasters + "\\cost" + str(uniqueID) + ".img")
        costRas = tempRasters + "\\cost" + str(uniqueID) + ".img"
        outCostDistance = CostDistance("uppoint", outRas, 500, tempRasters + "\\backlink" + str(uniqueID) + ".img")
        outCostDistance.save(tempRasters + "\\dist" + str(uniqueID) + ".img")
        distRas = tempRasters + "\\dist" + str(uniqueID) + ".img"
        backlinkRas = tempRasters + "\\backlink" + str(uniqueID) + ".img"
        outCostPath = CostPath("setpoint", distRas, backlinkRas)
        outCostPath.save(tempRasters + "\\path" + str(uniqueID) + ".img")
        pathRas = tempRasters + "\\path" + str(uniqueID) + ".img"

And here is the error:

Traceback (most recent call last):
  File "S:\Scripts\June6.py", line 138, in <module>
    outCostPath = CostPath("setpoint", distRas, backlinkRas)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 688, in CostPath
    destination_field)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 681, in wrapper
    destination_field)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 484, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 999999: Error executing function.
The table was not found. [CostPat_shp1]
The table was not found. [fras_bnd_CostPat_shp1]
The operation was attempted on an empty geometry.
ERROR 010029: Unable to create the raster C:\Users\amts\Documents\ArcGIS\Default.gdb\CostPat_shp1.
ERROR 010067: Error in executing grid expression.
Failed to execute (CostPath).

Both the backlink and cost distance rasters are saved and can be drawn in the map, so there doesn't appear to be anything wrong with them. Is this just 'one of those things' that happens in ArcGIS?


texts: Calculating area in polygons NOT covered by other polygons in ArcGIS Desktop?

texts:
I have two polygon shapefiles, one is county boundaries and the other is the National Flood Hazard Layer that shows areas in the floodplain. The flood layer is an extremely complex and challenging file to work with. It has many polygons that have 1 million + vertices and many overlapping polygons (that shouldn't be).

My goal is to find the area of counties either covered by the floodplain polygons or the area that is not covered by the polygons. Here's what I've done so far...

I first tried using identity to and then summary statistics on the county and flood layer. Unfortunately because there are so many overlapping polygons my output is wildly inaccurate sometimes with more than 100% of the county being covered by floodplain.

I then tried using topology with the polygons must not overlap tool and the create feature error fix. Unfortunately this crashes arcmap every time I run it.

My next thought was to calculate the area not covered by any polygon hence bypassing the problem of overlapping polygons. I'm currently trying to turn the polygon shapefile into a raster to do this but its a slow process.

I have computers running ArcMap 10.2.2 and 10.3 with advanced licences.


texts: How folded a widget on screen in Web App Builder?

texts:
I want make that a window like this

Has the foldable behaviour, like this foldable window:

I need a window foldable, moveable and resizable, any idea?


texts: Importing array directly as raster to PostGIS in memory?

texts:
Below is python code I wrote to import downloaded rasters into a postgis table using raster2pgsql. I have a question involving efficiency.

Instead of writing the raster to disk after doing array manipulations and then importing the written file to postgis, is there a way to directly import the array as a raster to postgis in memory? ie can I skip writing the raster to disk somehow?

I realize raster2pgsql is expecting a filename, but maybe there is a solution that goes straight from the in memory array to an in memory raster to postgis.

def write_raster(file_path, rast_array, no_data_value, rast_cols, rast_rows, projection, transform):
    driver = gdal.GetDriverByName('Gtiff')
    raster = driver.Create(file_path, rast_cols, rast_rows, 1, gdal.GDT_Float32)
    band = raster.GetRasterBand(1)
    band.SetNoDataValue(no_data_value)
    band.WriteArray(rast_array)
    raster.SetProjection(projection)
    raster.SetGeoTransform(transform)
    band.FlushCache()

def save_raster_to_postgis(raster_path, table_name):
    new_table = not table_exists(table_name)

    # remove old entry if already exists
    if not new_table:
        query = "DELETE FROM %(table)s WHERE filename = %(filename)s;"
        data = {"table": AsIs(table_name), "filename": os.path.basename(raster_path)}
        curs.execute(query, data)
        conn.commit()

    # insert the raster (either create a new table or append to previously created table)
    if new_table:
        import_command = "raster2pgsql -s 4269 -c -I -C -M -F -t auto {file} public.{table}"
            .format(file=raster_path, table=table_name)
    else:
        import_command = "raster2pgsql -s 4269 -a -M -F -t auto {file} public.{table}"
            .format(file=raster_path, table=table_name)
    import_command2 = "psql -h {host} -p {port} -d {database} -U {user} --no-password"
        .format(host=db["host"], port=db["port"], database=db["db"], user=db["user"])
    ps = subprocess.Popen(import_command, stdout=subprocess.PIPE, shell=True)
    subprocess.check_output(import_command2, stdin=ps.stdout, shell=True)
    ps.wait()

ds = gdal.Open(some_rast_file)
band1 = ds.GetRasterBand(1)
rast_array = band1.ReadAsArray()

# do some manipulations to the array here

write_raster(save_path + file_name + '.tif', temps_array, -9999, 2145, 1377, projection, transform)

save_raster_to_postgis(save_path + file_name + '.tif', 'my_rast_table')


texts: How to use ESRI 'Ecological Land Units' in QGIS?

texts:
I am trying to use the Map of Global Ecological Land Units –
An Ecophysiographic Stratification Approach (2014) in QGIS. I downloaded the World ELU 2015.zip from
https://rmgsc.cr.usgs.gov/outgoing/ecosystems/Global/

I am able to add the tif to the map and partially visualize the ELU with the default style classification QGIS (2.18.8) using grayscale pseudocolor.

ESRI has published a white paper with guidelines on on how to use the ELU on desktop but of course the explanation is ArcGIS and ArcGIS online based.  https://blogs.esri.com/esri/esri-insider/2016/04/06/updated-map-features-global-ecology-in-unprecedented-detail/

I never used raster data before.

So my question is on how can I:
1)obtain the information of dataset in QGIS (beyond just the image and pixel values by clicking on it), and
2)bring the ELU-ESRI color scheme to my QGIS map project


texts: Filtering stores location with login rules using GeoServer?

texts:
I want to filter the location/directory proposed by the stores, with the login rules.

More clearly, If my user "John" can only access/admin the "JohnWorkspace", I want he saw only one directory "dirForJohn" with the datastores for his Shapefile. and nothing else (no parent dir or what else).

Can I do that with GeoServer ?


texts: Using NLCD to quantify rate of urban sprawl in a watershed

texts:
I have downloaded five datasets and would like to quantify the rate of urban sprawl in a watershed from 2001 to 2011 using NLCD. Here is the data I have: (1) NLCD 2001 to 2011 Land Cover from to Change Index, (2) NLCD 2011 Land Cover (3) NLCD 2011 Percent Developed Imperviousness, (4) NLCD 2001 Land Cover (2011 Edition), and (5) NLCD 2001 Percent Developed Imperviousness (2011 Edition).

Once I clip the data to my watershed boundary, how do I calculate the rate of urban sprawl?

I know that I won’t need all of the above data, but figured this is essentially everything I need for my study.

All of the data was downloaded from https://www.mrlc.gov/index.php


texts: Write a PostgreSQL extension that makes use of PostGIS internal methods?

texts: As I've just found out, it is possible to extend PostgreSQL with custom functions written in C. And I also now know that that's how PostGIS itself is written. However what if I want to build on top of PostGIS methods, is it possible? So, like, write an extension that would make use of some of the PostGIS internal methods?


texts: Without Proxy page can't use IE8 and IE9 for Editing and Buffer

texts: The server where my editing application is hosted does not allow a proxy page because of outgoing calls restriction. This makes it difficult for IE8 and IE9 since using CORS is not an option as there is limited support for these browsers. The limitation is that they do not support POST when using XDomainRequest and also there is a URL length limit. What is the best practice to POST Cross Origin messages from these browsers. Even if I write custom code I would like to leverage the heavy lifting done by esri's tasks like the ApplyEdits event handlers.


texts:
I'm trying to upload a 20MB shapefile using the REST API, but gets a timeout error. I can upload smaller files without problems.

I'd like to be able to upload a shapefile of any size without restriction. If there's a way of specifying timeout, I guess it might be problematic as I wouldn't want to restrict file size; I do validation on the shapefile, and it should be uploaded provided it's valid.

Part of my code is shown below:

    private string getDbXml(string dsName)
    {
        string exposePKsParameter = "<entry key=\"Expose primary keys\">true</entry>";

        string xml = "<dataStore>" +
                        "<name>" + dsName + "</name>" +
                        "<enabled>true</enabled>" +
                        "<connectionParameters>" +
                            "<host>" + DB_HOST + "</host>" +
                            "<port>" + DB_PORT + "</port>" +
                            "<database>" + DB_NAME + "</database>" +
                            "<user>" + DB_USER + "</user>" +
                            "<password>" + DB_PASSWD + "</password>" +
                            "<dbtype>" + DB_TYPE + "</dbtype>" +
                            "<namespace>" + dsName + "_NS</namespace>" +
                            exposePKsParameter +
                         "</connectionParameters>" +
                      "</dataStore>";

        return xml;

    }

My environment setup: GeoServer 2.1.3, PostGIS 1.5, and C#.


texts: Slope Aspect Curvature algorithm giving error

texts:
I am trying to execute the slope aspect curvature algorithm in QGIS using a stand-alone script. When I run this code, I am getting a 'problem executing algorithm' error, which occurs when trying to run this slope aspect curvature algorithm.

I included my code for creating the clipped DEM and for running the slope aspect curvature algorithm.

processing.runandload('gdalogr:cliprasterbymasklayer',
                       raster_layer,    #INPUT <ParameterRaster>
                       mask_layer,      #MASK <ParameterVector>
                       "none",          #NO_DATA <ParameterString>
                       False,           #ALPHA_BAND <ParameterBoolean>
                       True,            #CROP_TO_CUTLINE <ParameterBoolean>
                       True,            #KEEP_RESOLUTION <ParameterBoolean>
                       5,               #RTYPE <ParameterSelection>
                       4,               #COMPRESS <ParameterSelection>
                       75,              #JPEGCOMPRESSION <ParameterNumber>
                       6,               #ZLEVEL <ParameterNumber>
                       1,               #PREDICTOR <ParameterNumber>
                       False,           #TILED <ParameterBoolean>
                       0,               #BIGTIFF <ParameterSelection>
                       False,           #TFW <ParameterBoolean>
                       "",              #EXTRA <ParameterString>
                       "C:/Users/Rachel/Desktop/clippedDEM.tif")     #OUTPUT <OutputRaster>


texts: Generate Statistics from check-out GDB While performing Check-in

texts:
I have created Check-out Replica from the main ArcSDE GDB. I have below questions.

How to generate statistics of the edits performed on the check-out GDB, before check-in. for example: i want to have a report on edits as compared to the parent version and location of edits, type of editd (Attribute, Spatial).


texts: QGIS - Rendering Order: Layer Order > Symbol Levels > Feature Rendering Order

texts:
A similar question was already asked two years ago, but no answers were provided, things may have changed, or if not developers may now consider this something for a future version.

Can I change the precedence when combining "Control feature rendering order" with QGIS symbol levels?

However that question was only about the last two affecting rendering order whereas my question also features data in multiple layers that also affects rendering order.

So it seems that there are three factors which possibly affect rendering order in QGIS and the order of preference is:

Layer Rendering Order
Symbol Levels
Feature Rendering Order

This doesn't seem to be officially documented and myself and the person asking the previous question have just discovered it while trying to get things to render in a certain order and finding it impossible.

The data I have is Ordnance Survey's Open Zoomstack. This breaks roads down into three separate layers: National, Regional and Local. Each of those layers also has a 'level' field to indicate bridges ('1') and bridges over bridges ('2') and so on, and naturally this is what feature rendering order is for.

Right away, the issue is that I can't seem to find a way to ensure that a level 1 feature is always drawn over a level 0 feature regardless of layer. I've ordered the layers National > Regional > Local and therefore National roads will always be drawn over Regional roads, even if the national road is level o and the regional road is level 1, meaning that the regional road should appear as a bridge over the national road.

Similarly symbol levels always override feature order (but not layer order). If you use two symbols to create a road with an outline (or perhaps three to create a road with a centre-line) and use symbol levels to ensure it looks right where roads intersect then firstly, you can only intersect roads in this way if they are in the same layer. So I can't have local roads intersecting with regional roads. Meanwhile where roads in the same layer should appear as bridges rather that intersections, then they will intersect anyway since the symbol level will override the feature rendering order.

I haven't used ArcGIS in quite some time but remember that you could make symbol levels apply to an entire group of layers in that. Perhaps this is a feature that QGIS should also have, as well as being able to do similar with feature order.

Perhaps a more flexible solution for QGIS would be if the render order for each individual symbol was a combination of symbol order, feature order and layer order. For example if the symbol order is 2, the feature order 3 and the layer order 1, then it will be drawn in rendering pass 6. If then you wanted to ensure that a certain feature was always drawn on top, you could set its feature order to something high like 100. If you wanted a certain symbol to always be drawn at the bottom then you could get the symbol order to -100 (or the expression in my case, 'level*100'), and so on. You'd need to be able to also assign numerical values to layer rendering order too.


texts:
I have trawled over many pages about the OpenLayers handlerOptions for drawing polygon holes. I can't create a polygon hole in existing polygons using the following:

var draw = new OpenLayers.Control.DrawFeature (
    farmpaddocks,
    OpenLayers.Handler.Polygon,
    {
      handlerOptions: {holeModifier: "altKey", multi: true, drawingHole:true,},
      title: "Draw Feature",
      displayClass: "olControlDrawFeaturePolygon",
      featureAdded: setAttributes
    }
);

All I get are tiny dots on the vertices... no polygons lines are drawn when pressing "alt" key.

Any ideas?


texts: Seeing all of image in Google Earth Engine

texts:
var s2SR_list = s2SR.toList(s2SR.size());

var s2SRMasked = s2SR_list.map(function (img) {

  return cld.sclMask(['cloud_low', 'cloud_medium', 'cloud_high', 'shadow'])(img);

});

print('s2SRMasked', s2SRMasked);

How can I see this 2SRMasked images?


texts:
I'm currently required to periodically pull cadastral lines from a WFS server and my current methodology is resulting in slow saves.
I first query a very large layer within the WFS with a known primary key to isolate 1 shape. I then adjust my map extent to cover the area I require and then Add the layer again except this time with no query and just the "Only request features overlapping the view extent" button checked.
This then shows me all the shapes within my view and I'm happy with that, the problem arises when I make a selection of say 20-50 polygons, right click the layer and then go export -> "Save selected features as.."
I choose DXF and then it hangs on 99%, the save does eventually work but it takes upwards of 15minutes and some random timeout errors. Once this save is done I can save another separate selection within the same screen instantly any ideas?
Also I would prefer lines instead of polygons so I usually do "polygons to lines" which also works but suffers the same 15 minute delay.

This is an example of the data and need in isolation instead of my entire state which will crash QGIS.

And this is where I'm trying to limit the geometry I get to my current map but with no success.


texts:
I am working in QGIS 1.8.0 Lisbosa. Currently I am trying to georeference a building blueprint in a .TIF format to a Google Streets layer. After selecting four points (the corners of the building), I run the function (Thin Plate Spline -  Cubic Spline) with a target SRS of UTM Zone 18N (EPSG: 3725). However, the results don't show up on the map. When I click zoom to layer, the screenview doesn't move. So either it isn't working, or the results are invisible.

Any idea what I'm doing wrong here? I know I've Georeferenced similar blueprints before, but I'm not sure what I'm doing differently.


texts: Adding polygon attributes (Shape_Area) to point data using ArcGIS ModelBuilder?

texts:
I converted my polygon feature class P1 stored in GDB to centroids by Feature to point tool. I would like to keep the Shape_Area and Shape_Length, native to original polygons, in my resulting Centroids.

I tried to add Shape_Area and Shape_Length to Centroid data by Spatial Join: the attribute table is copied, but the Shape_Area and Shape_Length are missing. I also tried the Intersect the data, without success.

Is there a way how to add attributes (Shape_Area, Shape_Length) from my original polygons to Centroids? I would like to keep the original polygon values the weight the centroids later.

I am working with ArcGIS 10.4, ideally in ModelBuilder.


texts: Reclassify a rasterList and sum the reclassified value from each Pixel

texts:
I'm trying to write a script in arcpy which reclassified a List of Rasters.
This rasters are made by another tool before so there is  already a loop. Actually I get the error 000582 and before this the script could reclassify but overwrote the rasters. In the end I like to count the pixel with the value 1 so that I create one raster with the sum of the value 1 as pixelvalue. How can I do this and what is wrong with the following script? Do I need the tool “plus” at the end?

import arcpy
import os
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
env.overwriteOutput = True
env.workspace = r"D:\Julia_T\projekt\aufgabe3\testdatenA"
Dir = env.workspace
#folder to list
ascFileList = arcpy.ListFiles("*.asc")
print ascFileList

for ascFile in ascFileList:
#geoprocessing steps
ascFileName = os.path.splitext (ascFile)[0]
print ascFileName
#works
rastFile = ascFileName + "_output.img"
print rastFile
#works
#run the tool
arcpy.ASCIIToRaster_conversion (env.workspace + "\\" + ascFile, env.workspace + "\\" +     rastFile, "INTEGER")
# Process: Projektion definieren
arcpy.DefineProjection_management (env.workspace + "\\" + rastFile,     "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRI    MEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
#list raster
rasters = arcpy.ListRasters("*", "IMG")
for ras in rasters:
    print "working..." + ras #works
    # new name
    name = "class_" + ras
    # Recalssify
    reclassField = "Value"
    remap = RemapValue ([["8",1],["2",0],["5",0]])
    outReclassify = Reclassify(ras, reclassField, remap, "NODATA")
    #save
    outReclassify.save (Dir + "\\" + name)

# Execute Plus
#outPlus = Plus(inRaster1, inRaster2)

I'm working with ArcMap 10.2.1.


texts: How can I convert esriMeters to Lat/Lng?

texts: I am pulling a JSON list of features into a web page. The points come as esriMeters and the rest of my app uses Lat/Lng. How can I do this conversion in the browser (javascript)?


texts: Geonetwork Authentication from an API client

texts:
I'm working on developing a client application that can push metadata to Geonetwork using Reactjs and Geonetowork API. To be more specific, I'm trying to add a record with the API /{portal}/api/0.1/records. My question is how to I can authenticate my request. I tried the Basic Authentication but I'm getting (403 Forbidden).
Sample of my request:
  let basicAuth = 'Basic ' + btoa("username" + ":" + "password");
  axios.put("/geonetwork/srv/api/0.1/records?metadataType=SUB_TEMPLATE", this.state.xml, {
    headers: {
      accept': 'application/json',
      'Content-Type': 'application/xml',
      'Authorization': + basicAuth
    },
  })


texts: Getting a GeoServer created style to Leaflet markers via PHP

texts:
How can I use a point/marker style created on GeoServer within my Leaflet map when I am drawing data from the server using PHP.
I am currently PHP to draw geodata from my GIS server, however I would like to use one of my custom styles created within GeoServer as the marker which shows on Leaflet as Leaflet does not support triangle or cross shaped markers by default. Is it possible to draw my marker type using WMS or PHP and apply it to the data drawn from my database using PHP?
Current marker options
var markerOptions = {
            radius: 4,
            fillColor: '#FF6600',
            color: '#000',
            weight: 1,
            opacity: 1,
            fillOpacity: 0.4 }

GeoServer custom style AKA Pub Style
*{
  mark: symbol(triangle);
  mark-size: 8px;
  :mark {
    fill:#FF0608;
    fill-opacity:0.4;
    stroke:#080;
  }
}

current call using AJAX to get PHP data using markerOptions as the style
var pTable = 'geoupload.pubs';
var pIRETable = 'geoupload.Pubs_and_Clubs_in_Causeway_Coast_and_Glens NI';
var pColumn =["id", "p_name", "address1", "postcode", "local_authority"];
var pGeom = "geom";

$(document).ready(function(){
    $.ajax({
        url: "knn.php",
        data:{table: pTable, inclu: pColumn, geomy: pGeom}
    })
    .done(function(data){
        pubsLay = L.geoJSON(JSON.parse(data),{
            pointToLayer: function (feature, latlng){
                return L.circleMarker(latlng,markerOptions);
        }});
        pub_map.addLayer(pubsLay);
    })
});

KNN.php call to get data from server
//call table and data
    $table = $_GET['table'];
    $inclu = $_GET['inclu'];
    $geomy = $_GET['geomy'];

    //construct substitutions
    $fields = "";
    foreach ($inclu as $i => $col) {
        $fields = $fields . $col . ", ";
    }
    $fields = $fields . $geomy;
    $geom = " - " . $geomy;

    //build query
    $sql1 =
    "SELECT jsonb_build_object(
            'type',     'FeatureCollection',
            'features', jsonb_agg(feature))
    FROM (
        SELECT jsonb_build_object(
            'type',       'Feature',
            'geometry',   ST_AsGeoJSON(st_transform(geom,4326),4)::jsonb,
            'properties', to_jsonb(row) #omits#
        ) AS feature
        FROM (SELECT #fields# FROM #table#) row) features;";

    //replace placeholders
    $sql2 = str_replace("#table#", $table, $sql1);
    $sql3 = str_replace("#fields#", $fields, $sql2);
    $sql4 = str_replace("#omits#", $geom, $sql3);
    ChromePhp::log($sql4);

//execute query
if (!$response = pg_query($conn, $sql4)) {
    ChromePhp::error("*SQL failed*");
    exit;
}

//return data
while ($row = pg_fetch_row($response)) {
    foreach ($row as $i => $attr){
        echo $attr;
    }
}


texts: How to add a perpendicular line to the start & end nodes of every polyline?

texts:
Is there a way to add a perpendicular line to the end node of each polyline?
Something like this:


texts: Georeferencing vector data to another vector in QGIS?

texts:
I have the latest version of QGIS Las Palmas 2.18 I am trying to georeference a vector data to another.

I saw most people use VectorBender plugin but I don't seem to find it in plugins


texts: Leaflet invalidate size - is not a function (undefined)

texts:
I have a leaflet map that is nested within a div.  I have tabs.  When the tab that contains the leaflet map is the default tab, the map loads and works correctly. When it is not, it does not load correctly.  It does not load all tiles.

If I resize the page, it then loads.  I am trying map.invalidateSize but get an is not a function error.  Leaflet object properties at the time the invalidateSize is invoked are as follows:

_leaflet_dblclickdblclick4_2:r()
_leaflet_events: {…}
_leaflet_id: 3
_leaflet_pointerdowndblclick4_2:n()
_leaflet_pointerupdblclick4_2:o()
_leaflet_touchstarttouchstart15_16()
_leaflet_touchstarttouchstart23_24()
accessKey: ""
accessKeyLabel: ""
align: ""
assignedSlot: null
attributes: NamedNodeMap(3)
baseURI: "http://localhost/RG/#10/44.1206/-122.3067"
childElementCount: 3
childNodes: NodeList(4)
children: HTMLCollection(3)
classList: DOMTokenList(6)
className: "leaflet-container leaflet-touch leaflet-fade-anim leaflet-grab leaflet-touch-drag leaflet-touch-zoom"
clientHeight: 345
clientLeft: 0
clientTop: 0
clientWidth: 1036
contentEditable: "inherit"
contextMenu: null
dataset: DOMStringMap(0)
dir: ""
draggable: false
firstChild: #text
firstElementChild: div.leaflet-pane.leaflet-map-pane
hidden: false
id: "map"

I have the following code:

function openPage(pageName, elmnt, color) {

 // Hide all elements with class="tabcontent" by default */
 var i, tabcontent, tablinks;
 tabcontent = document.getElementsByClassName("tabcontent");
 for (i = 0; i < tabcontent.length; i++) {
     tabcontent[i].style.display = "none";
 }

 // Remove the background color of all tablinks/buttons
 tablinks = document.getElementsByClassName("tablink");
 for (i = 0; i < tablinks.length; i++) {
     tablinks[i].style.backgroundColor = "";
 }

 // Show the specific tab content
document.getElementById(pageName).style.display = "block";

if ([pageName] == "Flow") {
document.getElementById("FlowLog").style.display = "block";
document.getElementById("lmap").style.visibility = "visible";
setTimeout(function(){ map.invalidateSize(true)}, 400);
}
}

 var map = L.map('lmap', {

//  preferCanvas: true,  //for leaflet image plugin
    measureControl:true,
    zoomControl:true, maxZoom:45, minZoom:1,
    center: ([44.120, -122.5951]),
    zoom: 11
   })


texts: metadata in Orfeo Toolbox Orthorectification

texts: I'm using Orfeo Toolbox to orthorectify a WorldView-3 Basic Imagery Product. The process seems very straight forward. I just input the image and a DEM of the area, and Orfeo gives me an orthorectified image, so apparently everything is very fine. However, I'm wondering whether the process is actually correct... how can I be sure that the software is actually using the information contained in the associated files (attitude file, metadata file, etc)?


texts: Map Grid Index / Reference for Canadian Topo Map

texts:
Does anyone know where I can find some KML at different resolutions to show the grid reference for the Canadian topographic maps (CanTopo and toporama). Much like is displayed on the map here: http://wms.ess-ws.nrcan.gc.ca/wms/indexcouverture_cantopo

As you zoom in the grids get more detailed so I imagine I'd require multiple KML files to provide the different levels of detail.


texts: XYZ tiles not rendering in QGIS 3.10 in some CRS

texts:
I recently installed the version 3.10 update of QGIS, and the xyz tile imagery stopped rendering. It originally showed up when opening a project saved under a previous version (3.8.3), but it also does not render if I add new layers.

I've tried adding several different xyz tiles (eg. Bing, Google, Open Street Maps, NAIP from web), and none of them show up. Local raster imagery shows up and renders fine. The xyz tiles render correctly if I open the projects in the older version.

I'm running QGIS(64-bit) 3.10.2-A Coruña on a Windows 8.1-pro 64 bit system.

After looking for similar problems on stack exchange, I tried running the script from this github link to connect with online resources following the accepted answer on question 285179. This added additional xyz layer options, but did not affect the behavior.  Layers are able to be added, but do not render under some projections/CRS's.
I tried starting with a clean, empty project and adding the xyz tiles before any local data, but this does not appear to affect the rendering issue.  If I add xyz layers in QGIS 3.10.2 and save the project, then open with a previous version of QGIS (3.8.3) they render fine.  This suggests it is not an issue with the layers not being added correctly.

After further investigation, the xyz tiles render correctly with some coordinate systems, but not others.
My original project was in EPSG 2285 (NAD83 Washington State Plane North). The tiles do not render in EPSG(ESRI):32148, EPSG 2286 (NAD 83 Washington South), EPSG 32048 or 32049 NAD27 Washington North or South, respectively.

These layers render fine in:
EPSG 4326, 3857 (versions of WGS 84)
EPSG 6339 (UTM zone 10N)
EPSG 6597 (NAD 83 (2011) WA state Plane)
EPSG 2855 (NAD 83 HARN Washington North)


texts: How to export only certain column from selected feature to csv file in QGIS using python code

texts:
I have selected features in QGIS, I want to export only two columns from the selected features to a csv file automatically using python code. I have tried the following code:

layer = qgis.utils.iface.activeLayer()
selected_features = layer.selectedFeatures()

for i in selected_features:
    idx = layer.dataProvider().fieldNameIndex('name')
    uv = layer.dataProvider().uniqueValues(idx)

for i in uv:
    if i:
        print str(i)


texts: Formula to convert State Plane Coordinates to Latitude and Longitude

texts:
The data I have is in State Plane Coordinates and I need them in regular Latitude and Longitude form. Is there an easy formula to convert them?

ZONE: 3104 (NAD 1983 StatePlane New York Long Island FIPS 3104 Feet)

I'm pulling the data from ArcGIS NYC MapPluto API query here
eg. Where: bbl = 1000430006 and I'm requesting GeoJSON output.

Is there a way to output the data already converted?


texts: Projecting shapefile with unknown CRS to WGS84 using QGIS?

texts:
I've tried setting QGIS to on the fly, changing the layer settings, nothing.

I'm at a lost cause. Although I am not a GIS professional.

http://www.ecaidata.org/dataset/spanish-and-mexican-land-grants-in-california/resource/9c7a68d8-a73e-40bf-8b39-a1f041735794?inner_span=True


texts: methods for solving this simple routing problem

texts:
In my transport geography class, we were asked to do a "best guess" on the route which does the following on this graph:

Starts at A
Visits each node
Returns to A
minimizes distance.

I wrote up program which calculated all of the five step paths that started from A, visited every node, and returned to A.  I then calculated the cost of these paths.

How would other people solve this?

My teacher gave us a "hint" that we should read about Dijkstra's algorithm, but I couldn't see how to apply that to solve this particular problem.  My solution would fall apart pretty quickly as the number of steps or nodes increased.


texts: What is the unit of measurement for scale?

texts: Maybe this is standard across mapping engines, but in MapGuide if I set the scale to 7 million I can see about half the world.  At about 4 million I can see a continent.  What unit of measurement could this scale possibly be in?


texts: I've been searching around to find NAIP data. The National Agriculture Imagery Program acquires aerial imagery during the agricultural growing seasons in the continental U.S. But I cannot seem to find much online for reliable data I recall there being NAIP data here: https://datagateway.nrcs.usda.gov/  but cannot seem to find it now.


texts: Adding information to polygons inside polygon from other polygon in QGIS?

texts:
I have two layers in QGIS. The first one contains buildings and the second one contains blocks. Now I want to add the ID of the block to all buildings inside this block.

Are there any functions in QGIS available for this?


texts: QGIS query syntax differs between attribute table and feature filter

texts:
Using QGIS 2.16.2, I can successfully select records from a shapefile attribute table using the following expression:

left("start_time", 10) = '2015-08-01'

Please note that start_time is a string field.

However, when I copy and paste the same expression into the shapefile's Layer Properties>General>feature filter>Query Builder I get the following error:

I would expect that the same expression would work everywhere in QGIS.  What am I doing wrong?


texts: GeoPDF with higher detail when zooming in in QGIS 3

texts:
In QGIS 3 there is a bult-in utility to create GeoPDFs but it has a limitation of 3000PPP max (resolution).

So when I open the GeoPDF everything seems fine:

But when I zoom in, the strokes are so wide that I can barely see anything....

Is there a way in QGIS 3 to build the GeoPDF so that when I zoom in I can see everything with a much higher detail?
Breaking the 3000PPP max. resolution limitation should be a solution I guess but I am not sure if it's possible...


texts: Per-pixel Descriptive Statistics through multiple Raster layers

texts:
I'm looking for computing different descriptive statistics through a set of raster layers. I need to make it with a python script capable of running into a standalone tool within the QGIS. To do so first I made it by using the "Model Builder" within the QGIS and from the model I have extracted the script. I used the r.series from GRASS and it works as a model. The script is the following, for the case of the "maximum":

##=name

##1=raster

##3=raster

##2=raster

##4=raster

##maximum=output raster

outputs_GRASSR.SERIES_1=processing.runalg('grass:r.series',[1,3,2,4],True,6,'-10000000000,10000000000',None,0.0,maximum)

Actually at this point I need to use a "multiple raster" instead of 1,2,3,4. I mean I would like to make something like:

##MAXIMUM=name

##RASTERS=multiple raster

##maximum=output raster

import os
import processing

lddLrs = RASTERS.split(';')
for lyr in lddLrs:
    processing.runalg('grass:r.series', RASTERS,True,6,'-10000000000,10000000000',None,0.0,maximum)

But this doesn't work. Any idea to solve it?

I have tried the following:

##=name

##1=multiple raster

##maximum=output raster

import processing

processing.runalg('grass:r.series', [1],True,6,None,0.0,'-10000000000,10000000000',maximum)

But it still doesn.t work. Any solution?


texts: Display coordinates of points in DMS format in print composer

texts:
In QGIS 2.14.2 I have a series of polylines of which I want the start and endpoint to be labeled with its coordinates in WGS84 (EPSG 4326) like this: 51*15'23.40"N, 5*17'56.58"E

I thought I figured a workaround by:

Creating new point shapefile and creating points manually at start and endpoint of each polyline.
Open field calculator and create a virtual field "X" containing $x and another containing $Y (Decimals)
Label them by X||'\n'||Y

However the field calculator gives me an output like this
X:5.26344964931 Y:51.20551633106

How do I display them in DMS format?


texts:
I have some lines which need to be sliced by a certain portion.
I can get the outcome by using turf.lineSliceAlong function, something like:
totalLength = turf.length(myLine);
markerSliced = totalLength * (1 - leftPortion);
slicedLine = turf.lineSliceAlong(myLine, markerSliced, totalLength);

I'm wondering is there a better way of doing this by using pure Leaflet functions?


texts: How to construct semi-circle-shaped areas of confidence in ArcGIS Desktop?

texts: We have polygon maps of river-bottom type and are interested in determining if particular bottom-feeding fish prefer one bottom type over another.  We locate an acoustically tagged fish at a point in time using a boat with an acoustic array. We have a measure of confidence of the location, but the measure of confidence is not a circle; it is more pie-shaped; i.e. we know that the fish is located in front of the boat in a semi-circle of know radius.  We would like to create a semi-circular polygon at each location point with the semi-circular diameter perpendicular to the boat heading.  We would like to use the resulting semi-circular polygons to clip out pieces of the underlying bottom type map. I am looking for a procedure in ArcGIS that would allow me to create these oriented semicircles given a table of locations and boat headings.


texts:
I am producing some maps of a point layer using arcpy.mapping.ExportToPDF. Is there a possibility of creating HTML popup for the final PDF based on this link?

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Using_HTML_pop_ups_for_feature_layers/00s500000018000000/


texts:
I have installed GDAL using conda and it installed successfully. But when I try

import gdal in jupyter notebook, it is showing this error.

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-25-ee38efafc30b> in <module>
----> 1 import gdal

~\Anaconda3\lib\site-packages\gdal.py in <module>
      1 # import osgeo.gdal as a convenience
----> 2 from osgeo.gdal import deprecation_warn
      3 deprecation_warn('gdal')
      4
      5 from osgeo.gdal import *

ModuleNotFoundError: No module named 'osgeo'

and when I try from osgeo import gdal, it is showing this error:

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-26-61e0eefce4de> in <module>
----> 1 from osgeo import gdal

ModuleNotFoundError: No module named 'osgeo'

for GDAL installation i used this command:

conda install -c conda-forge gdal

version information

-- python 3.7
-- conda 4.7.5
-- gdal 2.3.3

How to solve this?


texts:
I have a corel draw (.cdr) file that I would like to import into a GIS. It is a vector format which I thing has no geographic information associated. I would like to be able to get it into MapInfo Professional and use the vector registering tools to place it correctly.

The problem I have is I cannot get it into .tab easily.

Does anyone know a way of getting it into any GIS/CAD format I could work with

Thanks,


texts: Getting RasterSymbolizer from a StyledLayerDescriptor?

texts:
We are using Geotools for a Java GIS project. I have an org.geotools.styling.StyledLayerDescriptor object created from a XML stream and would like to create an org.geotools.styling.RasterSymbolizer from it.

How do I go about to do this?

The documentation (http://docs.geotools.org/latest/userguide/library/render/index.html) does not give any examples for raster symbolizers.


texts: Non-square raster pixels created during TIN Interpolation?

texts:
I am trying to create a raster by using the TIN Interpolation tool in QGIS. I'm using a 3D polygon shapefile as input. I'm looking to create a raster of square pixels, although despite specifying a square output extent and 0.1m cell size, the resulting tif file contains cell size that are slightly non-square (e.g. 0.09999500000000000055,-0.09999750000000000305).

I've done this process many times before with no issue and can't understand why this is different. Can anyone tell me why this is happening?


texts: is there an operation for measuring constancy (ecology)

texts:
I'm looking for an effective way to measure changes in spatial location in a time series data set (whether an observation repeats in the same location between successive time-steps, or if it shifts and by how much).

What I have is a point file, showing the presence/absence of a species in four stages, in study region A.  I have a similar dataset for four other regions (B-D). I am looking for an approach to measure how much movement there is in region A over the four time-steps, and if this is greater, lesser, or equal to comparable movements in regions B-D.  Each region can be run as a separate analysis.

I am using ArcGIS (10.1 and 10.2), as well as GME.  I don't find any out of the box solution to what I'm trying to show.  The closest is to use a summary measure of point-distance.  The closest theoretical approach I have found is Colwell's (1974) concept of temporal/spatial constancy, used in Ecology.

If anyone has suggestions or knows of a tool in ecological modeling for measuring constancy, this would be very helpful.


texts: What's the mature solution for label placement?

texts:
I am working on an android map application, while I need to place the labels for points/lines/polygons manually.

I have googled and I found this:

https://www.ads.tuwien.ac.at/research/labeling/

But it seems that it is last updated at 2004, I am not sure if this manner is out of date. And further more, it just lable the points, while I need to other kind of feature like line,polygon.

So I wonder if there is a better solution?

Maybe I do not make myself clear enough.

I am not looking for what should be placed in the map but how to make the features in the map labeled easy readable and without overlap.


texts:
I'm trying to make a little tool who do a concave hull, and calculate its surface in a new field.
I'm not a python pro, so I began with treatment modeler like this :

It works very well !
But I have to treat more than hundred shapes, so I decide to save my model in python script and edit it to treat all shapes in a folder.
In first time, I wanted to try python script I saved. But strangely, it doesn't work... I have a weird message...

Here is my code :

##Calcule de surface cultivee=name
##couchedepiedsdevigne=vector
##couchedesortie=output vector

outputs_QGISCONCAVEHULL_1=processing.runalg('qgis:concavehull', couchedepiedsdevigne,0.3,False,False,None)
outputs_QGISFIELDCALCULATOR_1=processing.runalg('qgis:advancedpythonfieldcalculator', outputs_QGISCONCAVEHULL_1['OUTPUT'],"SURFACE",1,10.0,3.0,True,'$area',couchedesortie)

I don't understand why it's working in treatment modeler and not in python script...
Somebody can help me ?


texts: How can I assign an id column in SQL query that merges geometries?

texts:
I have a layer that consists of a number of features.  Each features is a line segment and has an attribute named zone_name which takes one of 7 values.  What I want to do is to use the atlas generation to print out a map for each zone (which consists of multiple features/lines) but (assuming I understand things correctly) the atlas generator iterates over all features in the layer -- first time I tried I got a large number of very small maps ;)

After a bit of research I came up with this:

SELECT zone_name,
     ST_Multi(ST_Collect(f.the_geom)) as singlegeom
     FROM (SELECT layer_name, (ST_Dump(the_geom)).geom As the_geom
                FROM
                tracks_master, layers where feature_gid = gid ) As f
GROUP BY zone_name

in the database manager SQL window.  When I run the query it produces a table with a single feature for each value of zone_name.  The only problem is that there is no integer id field so I can not load this as a layer.  I.e. when I select "load as layer" there isn't anything valid to put in "column with unique integer values".

There must be a simple way of doing this but my SQL isn't up to it.


texts: What are some alternative sources for USGS data?

texts:
I've been using the USGS National Map viewer (http://viewer.nationalmap.gov/viewer/) for getting data for QGIS, but the download system is driving me insane. I enter my email address and [sometimes] get a notice that my data will be available in up to 2 days (why so long???) If that email does arrive, occasionally the URLs are linked to a 0 byte file and I'm out of luck.

Is there any other place to get the data? This "click and hope your data arrives in 2 days" system keeps getting me in trouble and behind on my work. I don't understand why I can download multiple GBs of elevation data from USGS on demand but hydro data that's a fraction of the size requires a waiting period.

Edit: I'm specifically looking for hydro, transport and land cover for Kittitas Valley in Washington state. Thanks for the great answers.


texts: QGIS multiply(add) two raster but divide(deduct) another raster

texts:
I am making a parking analysis in one town. I need to obatin the zones where the parking of car is harder because there is less places to park and higher demand of cars.I have three rasters.

First raster represents the number of PUBLIC car parking places in a town (the redder the more parking places):

Second raster represents the number of PRIVATE car parking places in a town (the redder the more parking places):

I need to add (multiply) these two layers. To do so, in blending mode i select Multiply in both layers, and i get:

Which looks great. I'll call the union of these two layers A from now on.

The problem comes when a third layer comes into action, B from now on.
This third layer represents the vehicles in the census of the town, which is, the demand of car parking places:

What i think i need to achieve is B-A so that the places with high park demand (many vehicles in the census) in B, are lightened by A.

I don't know if I explained myself.

Regards,


texts: Replacing multiple values in attribute table field using ArcGIS field calculator and python parser?

texts:
I'm new to python and trying to write a python script using the field calculator in ArcMap for a specific field in an attribute table that will replace some of the values in the field but leave the other values as they are within the table.

For example: ""fourth" must be changed to "4th" whereas "Neilston" should remain as it is in the field.  The field in the table looks like this:

The python script used in the codeblock looks like this:

def streetNUM(field):

   if field == "fourth":

      field.replace ("fourth","4th")

   elif field == "fifth":

      field.replace ("fifth","5th")
   else:
      return field

And then in the box below the codeblock I have:

streetNUM(!FULLNAME!)

Nothing seems to change and I get an error code saying:

  Field is not nullable [FULLNAME]

Any suggestions as to where I may be going wrong with my code?


texts: Geonetwork 4 - Error accesing Swagger API

texts:
I recently installed Geonetwork 4.0.0.0 on my server, and when I try to access the Swagger API (http://XXX.XX.XXX.XX.:8080/geonetwork/srv/api/doc)

.
I am getting this error:

What is missing on my installation for access the API?


texts: What controls the direction of TopoGeometry Linestrings?

texts:
In PostGIS, what if anything controls the direction that topogeometry LINESTRINGs will have when converted to geometry?

It doesn't appear to be the order in which the edges are listed in the call to CreateTopoGeom. Which makes sense in the presence of loop edges, I guess, although presumably you could define the direction of a linestring that starts with a loop edge to be in the direction of that edge.

Is it even stable? If I do a bunch of edge splitting/healing, is it possible that the system will flip my linestrings that include those edges? If it's at least guaranteed to be stable I could include a should_be_reversed boolean column in the table that includes the linestring topogeometry column, and call ST_Reverse on the resulting geometry in a view?


texts:
I set up a Web Processing Service with 52°North and GRASS (using the wps-grass-bridge). All went smooth and I get a proper DescribeProcess for i.cluster. As you can see, I need to give the parameter "classes":

I installed Horst Duester's WPS plugin for QGIS and wanted to test my WPS now. Surprisingly I get all paramters as fields to fill beside of the obligatory "classes". It is just missing.
Does anybody have an idea what is wrong here?


texts: What causes Scripting Error 000714 in ArcGIS when using GetParametersAsText?

texts:
I am trying to create a script tool in ArcGIS 10.2 in order to incorporate into a model. The script needs to incorporate three inputs which i've built using GetParameterAsText function then setting up the corresponding links within ArcGIS Script Tool Wizard.

However, when I run the script I get ERROR 000714 (see below). Not sure how to handle this error or why it isn't working to begin with.

import arcpy

inFC = arcpy.GetParameterAsText(0)
outFCName = arcpy.GetParameterAsText(1)
outFolder = arcpy.GetParameterAsText(2)
outFC = outFolder + "\\" + outFCName

try:
    if arcpy.Exists(outFC):
        arcpy.Delete_management(outFC)
    arcpy.CreateFeatureclass_management(outFolder,outFCName,"POLYLINE","#","DISABLED","DISABLED",inFC)
    arcpy.AddField_management(outFC,"inFID","LONG","#","#","#","#","NULLABLE","NON_REQUIRED","#")

    iCursor = arcpy.da.InsertCursor(outFC, ["inFID","SHAPE@"])

    print "Starting!"

    with arcpy.da.SearchCursor(inFC,["OID@", "SHAPE@"]) as sCursor:
        for row in sCursor:
            inFID = row[0]
            # Print the current multipoint's ID
            #
            # print("Feature {0}:".format(row[0]))
            partnum = 0

            # Step through each part of the feature
            #
            for part in row[1]:
                # Print the part number
                #
                # print("Part {0}:".format(partnum))

                # Step through each vertex in the feature
                #
                prevX = None
                prevY = None
                for pnt in part:
                    if pnt:
                        # Print x,y coordinates of current point
                        #
                        # print("{0}, {1}".format(pnt.X, pnt.Y))
                        if prevX:
                            array = arcpy.Array([arcpy.Point(prevX, prevY),
                                                 arcpy.Point(pnt.X, pnt.Y)])
                            polyline = arcpy.Polyline(array)
                            iCursor.insertRow([inFID,polyline])
                        prevX = pnt.X
                        prevY = pnt.Y
                    else:
                        # If pnt is None, this represents an interior ring
                        #
                        print("Interior Ring:")
                partnum += 1

    del iCursor

except Exception as e:
    print e.message
    arcpy.AddError(e.message)

print "Complete"

Here is the error:

Executing: SplitLines "F:\Extra Data\Prop Frontage\test_for_py.shp" tester.shp "F:\Extra Data\Prop Frontage"
Start Time: Thu Jan 23 11:38:46 2014
Running script SplitLines...
ERROR 000714: Error in script SplitLines.
Error in executing: cmd.exe /C F:\MASTER~2  "F:\Extra Data\Prop Frontage\test_for_py.shp" "tester.shp" "F:\Extra Data\Prop Frontage"

Failed to execute (SplitLines).
Failed at Thu Jan 23 11:38:46 2014 (Elapsed Time: 0.16 seconds)


texts: Performing Raster Noise Reduction and Edge Smoothing?

texts:
I have a raster image like this

I want to delete the isolated noise point and smooth the edge (say the white single points, and I want the edges of the entity look smooth). How can I do this in ArcGIS or in R?

The cell size is 30*30.


texts: Viewing or importing *.pbf file (OSM data, binary format) in ArcGIS Desktop or QGIS?

texts:
How to view or import a pbf file (OSM data,binary format) in ArcGIS Desktop or QGIS?

Are any 3rd party applications needed?

I would like to edit this file on any GIS Software.


texts:
I'm using QGIS and I have some large features that need to be diced up into features with less vertices.

If I were using arcGIS I could use the dice tool

If I were only interested in lines I could incorporate mapplus' solution here

Is there an equivilent tool that I could script in pyQGIS?


texts: Using flags in new-style QGIS Processing Algorithms

texts:
In QGIS 3.6+, there are two different ways to construct a Processing Algorithm.

Defining a class inherited from QgsProcessingAlgorithm. In this case, I can define flags for an algorithm as follows:
class ExampleProcessingAlgorithm(QgsProcessingAlgorithm):
    ....

    def flags(self):
        return super().flags() | QgsProcessingAlgorithm.FlagNoThreading
    ....

Using @alg decorator.
@alg(name='algorithm_name', ...
@alg.input(type=alg.SOURCE, ...
@alg.output(type=alg.NUMBER, ...

def algorithm_name(instance, parameters, context, feedback, inputs):
    """
    Description of the algorithm.
    """

    ...

In @alg case, how can I set flags? There is no explanation about this in QGIS Documentation.


texts:
I'm working in Mapbox Studio (Classic) and want to prevent labels from rendering more than once.

I'm using text-repeat-distance to prevent labels from showing more than once:

#Admin_level3_Boundary {
  line-width: 1;
  line-color: #aaa;

  [zoom>=8] {
    text-name: [AREANAME];
    text-face-name: @sans;
    text-fill: #ccc;
    text-size: 8;
    text-avoid-edges: true;
    text-repeat-distance: 50000;
  }
}

According to the CartoCSS attributes outlined here, a high text-repeat-distance should do the trick.  But it's not working:

I've tried using text-min-distance and text-margin as well, but it causes other (non-repeat) labels to disappear if they are close by.

How can I prevent repeat labels from showing?


texts: How do I access GeoTransform array from gdal on the command line?

texts: I'm using the gdal library via OSGeo for windows, and I was just wondering how to access the geotransform array for a .ecw file.


texts: Show/Hide labels of Vector Tile layers on button click

texts:
I have displayed vector tiles layers with "application/vnd.mapbox-vector-tile" format in Open Layers 3.While creating vector tile I have applied only image style. I want to display labels for existing vector tile layer when user click on button.

While creating vector tile layer I am able to display labels because we will get features so I can apply text style. But in my case when user wants see labels they click on button that time labels has to display.

How to achieve that?


texts:
I work with many different projections and often times find my self working with a projection I do not recognize.

For example, I often times will get a generated CRS I do not recognize

(* Generated CRS (+proj=lcc +lat_1=37.08333333333334 +lat_2=38.66666666666666 +lat_0=36.33333333333334 +lon_0=-85.75 +x_0=1500000 +y_0=999999.9999999999 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs))*

Instead of searching for this projection I would like to import the projection from the layer. Arcmap does this in the project tool where you are able to import projections from a layer.

Is there a way to do this in QGIS?


texts:
I have a layer I need to symbolise with graduated colours based on one field and graduated symbols based on another field.

Sizes are based on a field that stores a count (range 1 to 50).  Colours are based on a field that stores cost (range 0 to 7500).

I set up Symbology by selecting Multiple Attributes, choosing my Count field and manually setting the sizes, then selecting Variation by Colour Ramp and choosing my Cost field for the ramp.

It has mostly worked, except it's not using quite the same colour ramp as I selected.  The two highest values are using the same colour, and I can't find a way to update it.

Symbology > Multiple Attributes:

Colour Ramp:

End result (note the high values are still yellow, not red):

The largest point on the map in my screenshot above should be coloured red, as the value for this point is over $4000.

Has anyone had success in setting symbology this way?  What do I need to do/change to get this to work?


texts: QGIS Sum Line Length Problem: Values 'nan'. Road length in Grid

texts:
I wanted to use the QGIS function 'Sum Line Length' to compute the sum of the road lengths in each grid cell in Nigeria.
My roads data is from here: https://data.humdata.org/dataset/nigeria-roads
My grid net is the PRIO Grid, hence 0.5 x 0.5 degree grid cells.
When I first used the function mentioned above all the data was in the "standard" geographical CRS EPSG:4326 - WGS 84. So the units were degrees. But since I need meters for better interpretation in my (spatial) regression, I exported the data to a projected CBS (EPSG:3786 - World Equidistant Cylindrical (Sphere)) and reran the function. Before when everything was in degrees, the column 'length' had actual values/numbers. But when I did it in the other CRS every cell which should have had values, has 'nan'.
I already tried different road data and different CRSs, nothing helped.
Although a similar question has been asked (QGIS "Sum Line Lengths" return 'nan.'), I don't see the solution.
(using QGIS 3.12.2)


texts:
I have a GeoJSON which returns a set of points, for each point I have different azimuth(rotation) is required.
I have tried using the rotate property of OpenLayers, however when I use it all the point in the GeoJSON are getting rotated to the angle defined.
t.setStyle(
      new ol.style.Style({
          image: new ol.style.Icon({
              src: "/Icons/signal.png",
              anchor: [ 0, 0.5 ],
              rotation: 0,
              scale: [ 0.35, 0.35 ]
          }),
      })
   );

Instead, I want to rotate each point in GeoJSON at a different angle by the azimuth value which I have.
Could someone suggest me a solution for this?


texts: Using Field Calculator to assign consecutive numbers to a field based on attributes in another field

texts:
I am trying to assign consecutive numbers to a numbers field using Field Calculator so that all my 500 features have unique label values, however some features have multiple polygons and therefore have two or more rows in the attribute table. in this case I need those features to have the same numerical value. i have another field with the names of all features that i would like to base this numbering system off of. So, i need a unique label for all unique features starting from 1 to the end.

How should i go about doing this?


texts:
A picture is worth a thousand words:

(source: a traffic modelling suite) This is what I want to create in Qgis: a plot where the thickness of a line depends on the sum of three attributes (let's call them A, B and C) where the bandwidth of each attribute depends on its value and has its own colour.

I know how to make the line thickness variable depending on the value of an attribute, but is it possible in Qgis to combine different attributes like above?


texts:
I'm trying to add a GeoJSON polygon layer from a file, give it a specific style and show the label (ToolTip) but couldn't find a way to do it in a single function call
.
I only found different code samples doing each task alone (styling data or showing labels), then I mixed both codes and got something weird, it just created another layer on top of the previous one.

  var geojson = L.geoJSON(myGeoJson);

  L.geoJson(myGeoJson,{
    onEachFeature: function (feature, layer) {
        layer.bindTooltip(feature.properties.CASOS.toString(), {permanent: true}).openTooltip();
    }
    }).addTo(map)

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

How could I show labels and style the polygons without creating a stack of layers?

EDIT: I edited the question because I think I was unclear, I want to style the data, not the labels.


texts:
I have a raster file and I want to perform the method "polygonize". Then the problem occurs, that is shown on the screenshot.
It doesn t have anything to do with the rasterfile format, I veryfied.
Other functions in the same menue(Raster->Conversion do not work either, because of the same problem).
I m using the newest version of Qgis (1.8).
Older threads here couldn t (as far as I saw) help either.
So: I want to ask if anybody has an idea how to solve this problem?
Thanks.


texts:
I want to have a gridPanel with feature geometry. But I can't image how to do this. I have a FeatureStore

st = new GeoExt.data.FeatureStore({
    layer: myVecLayer,
    fields: [
        {name: 'id', type: 'int'},
        {name: 'mydata', type: 'String'}
    ],
    proxy: new GeoExt.data.ProtocolProxy({
        protocol: new OpenLayers.Protocol.WFS({
            version: "1.0.0",
            srsName:"EPSG:900913",
            url: "http://localhost:8080/geoserver/wfs",
            featureType: "filedata_temp",
            featureNS: "http://www.opengeospatial.net/cite",
            filter:Ffilter
        })
    }),
    autoLoad: true
});

How do you get the feature geometry (with separate X and Y) from this?

I know that I'm going to get WKT from FeatureStore and parse it to get the coordinates, but what next, use another store to hold the coordinates?


texts: Use MeteoFrance WMS with WebApp Builder as a WMSLayer load by a widget

texts:
I try to use WMS OGC services for weather forecasting I'm able to load WMS layer but when I move on the map layer follow the moving as an image who is stuck on the top of my screen. I set all the parameters and I can't figure out where is my problem...

Here is my code:

_loadMeteoFranceData : function (MFApiKey){
var requestMapOpts = {
    version : "1.3.0",
    LAYERS : "TOTAL_WATER_PRECIPITATION__GROUND_OR_WATER_SURFACE",
    CRS : "EPSG:4326",
    styles : "EAU__GROUND__RADAR_SHADING",
    format : "png",
    width : "2048",
    height : "2048",
    BBOX : "38,-12,53,8",
    dim_reference_time : "2017-07-04T03:00:00Z",
    time : "2017-07-04T04:00:00Z",
    token : MFApiKey
};

var urlGetMapArome = "https://geoservices.meteofrance.fr/services/MF-NWP-HIGHRES-AROME-0025-FRANCE-WMS";
var WMStest2 = new WMSLayer(urlGetMapArome,{
            format : requestMapOpts.format,
            version : requestMapOpts.version,
            visibleLayers : [requestMapOpts.LAYERS],
            opacity : 0.85,
            visible : true,
            resourceInfo: {
                getFeatureInfoURL: urlGetMapArome,
                getMapURL: urlGetMapArome,
                layerInfos: [
                      new WMSLayerInfo({
                name: requestMapOpts.LAYERS,
                title: "Quantité de précipitations sous forme liquide.",
                spatialReferences : new SpatialReference(4326),
                allExtents : [new Extent(38,-12,53,8, {wkid: 4326}),
                     new Extent(4230140.65014,-1345708.40841,5899933.01204,893463.751013,{wkid: 3857})],
                extent : new Extent(38,-12,53,8, {wkid: 4326}),
                legendURL : "https://geoservices.meteofrance.fr/api/"+MFApiKey+"/MF-NWP-HIGHRES-AROME-0025-FRANCE-WMS?service=WMS&version=1.3.0&sld_version=1.1.0&request=GetLegendGraphic&layer=TOTAL_WATER_PRECIPITATION__GROUND_OR_WATER_SURFACE&style=EAU__GROUND__RADAR_SHADING&format=image/png",
                description : "QTé précipitations modele AROME"
                      })],
                extent : new Extent(38,-12,53,8, {wkid: 4326}),

                customLayerParameters : {
            styles : requestMapOpts.styles,
            BBOX :requestMapOpts.BBOX,
            CRS : requestMapOpts.CRS,
            dim_reference_time : requestMapOpts.dim_reference_time,
            time : requestMapOpts.time,
            token: MFApiKey
            }}
            });

_map.addLayer(WMStest2);


texts: How stable is MapWindow on a 64 bit machine

texts:
One of my favorite ArcGIS extensions (N_SPECT) is no longer going to be supportd under ArcGIS 10. It will be totally rewritten as a MapWindow plug-in, so I need to get up-to-speed with that in case a project comes along where I need to use the tool.

I have a 64 bit Windows 7 machine, but the MapWindow website lists the 64 bit version as Experimental.
Has anyone here used the 64 bit version?
Have you found it to be less stable than the 32 bit version?

To clarify - I'm looking at the desktop application and use of associated plug-ins. I'm an end-user, and will not be working on or developing anything with the MapWindow Active X Control. The N-SPECT tool was process intensive as an ArcGIS extension, and I assume it will be in MapWindow, too. That's why I'm looking at the 64 bit version. I suppose all of this may be moot depending on how NOAA develops their new N-SPECT tool?


texts: I need Sentinel 2 imagery of Nigeria, I signed up in Copernicus Scihub but typing Nigeria into the text box no results appear. I must be doing something wrong, how can be possible that they do not have images for Nigeria?


texts: LecoS batch processing - Is there an easier way?

texts:
I'm using the LecoS software: https://github.com/Martin-Jung/LecoS

I've clipped my raster to reduce processing time to focus in on the areas I'm interested in. The entire raster has 11 classes, but if I clip the raster some of the classes also get clipped. I learned to look at the histogram of the different classes for each clip. If a class is missing, then I select "no" under the classification settings. My reason for retaining the entry in my batch option is that there is no way to delete individual rows - deleting removes rows from the base. However, there are some cases where the software runs into issues for some classes using particular metrics. For example, Euclidean Nearest-Neighbor Distance runs into problems for some classes and this has to do with the way that it is calculated, but there is no way to know this from the onset.

My batch processing involves 144 line items. It is painful. The batch processing window has an option to save the entries, but it defaults back to the drop-down menu options. For example, if a processing error occurs I have to re-open and it defaults back to Raster Classified: "Yes", Metrics (single class): "Landcover", and Load in GIS: "Yes". This means that I have click back through each drop-down item 144X for each column to align with options. It took me a long-time to individually enter the baseline options, but it was faster than doing the analysis over and over again. It is a process of discovery in the analysis on what items will not compute (e.g., a Euclidean metric for Agriculture 500 m buffer overlay - raster cell class 8 freezes in one of my smaller clips, but not others and the 1000 m buffer overlay works just fine). I have to note where it freezes - select that item out - then start over again.

I am wondering if there is an easier way to do this - I'm sure there is. I was hoping that someone might be able to direct me to a means to enter into the "back-door" using Python code. I don't know Python very well - but I'm reasonable with r-stats and a few other languages. This batch process is so clunky - it also saves each metric into separate files that I have to piece back together using r-stats script. It would also be nice if I could run a complete analysis including max, min, sd, eveness, etc. of the metrics, but this adds yet another layer of complexity in clicking dropdown menus and having to later re-integrate the separate csv files. Hoping someone could direct me to a scripting option for this type of analysis. I'm willing to learn the Python basics, but hoping someone could get me started. Is this a common issue with batch processing analysis in QGIS?


texts: Zonal statistics projected differently ArcGIS Pro 2.5

texts:
I have a fishnet of equally sized raster cells of 1/4 degree, which spans the whole world. Furthermore, I have a netCDF file containing precipitation data. I load both, the shape and the raster file, into ArcGIS Pro 2.5 and both have the same projection (WGS 1984 EPSG:4326). When I run Zonal statistics (ignoring no-data fields), I get the following result:

Grey areas are those without precipitation. As one can see, the newly created zonal statistics layer is shifted downward with respect to the original precipitation layer (the blueish map below), however the field entries seem to be correct (Areas with a lot of precipitation like Greenland, Brazil or West Canada are also shown as areas with high precipitation in the zonal statistics layer). The following picture shows the original precipitation layer to more easily compare the two.

My question is: Why is the zonal statistics layer shifted downward? Should I be concerned about the correctness of the zonal statistics being calculated? How, if necessary, can I fix this problem? (Note that the fishnet has another extension than the precipitation data. The antarctic region e.g. is not part of the fishnet)


texts: Adding latitude and longitude around an image saved in ArcMap

texts: Is there a way to add latitude and longitude around an image in ArcMap when you are in layout view? I have attached an example of what I mean. I would like to do this so people looking at a saved image of a map can know the general location of the image.


texts: Which features does QGIS team test when shifting from older to newer versions?

texts:
Which features in QGIS are covered by automatic tests of the QGIS team?

Ideally you could say most functions (such as add vector layer, saving a selected vector feature via save as selection) don't need user testing between versions and should just work.  However new features are added and things/processes are changed between versions.

This - I assume - would leave any custom workflow processes to be tested, particularly in commercial environments. And users should be using long term release (LTR) QGIS versions.


texts: Reducing arrays within arrays using Google Earth Engine?

texts:
I have a very large array containing lists, each with 364 values of daily cumulative rainfall for the years 1988-2018. I would like to calculate a mean for each day in this array, with the output being 364 values in an array, each value containing the mean for that day. Below is an example of the structure of my data and my desired outcome:

var array1 = [1,3,6,8,12,24,4,2,10,11,13]
var array2 = [5,9,2,18,10,2,6,8,4,3,3]

var arrays = [array1, array2]

var arrayMean = [3,6,4,13,11,13,5,5,7,7,8]

The arrayMean variable contains a list of values equally long to the input arrays, but with the mean for each value.

How do I perform this type of calculation?


texts:
I am using QGIS 3.6.3 and I am an amateur QGIS user.
My QGIS project currently has a map outlining suburb boundaries in New South Wales, Australia. I am wondering if I can create a buffer around particular suburb boundaries? If yes, how would I go about this?


texts:
I use SQL Server (2012 sp3) DB with Geoserver and i noticed that Geoserver likes to call SQL like

SELECT ... geom.STAsBinary() FROM (geoserver sql view) WHERE geom.filter(...)

In my test case .STAsBinary() triples query times, so i was wondering that is there way to to just read WKB from database

Like

 SELECT ... geom_wkb from (select ... ,geom_wkb, geom ) where geom.filter(...)


texts: I got the tiles of a raster and I would like to create the original raster

texts:
I got the tiles of a map named

txt
2x0x0.png
2x0x1.png
...
2x1x1.png
2x1x2.png
...
2x45x1.png
2x45x2.png

The first part is the zoom, the second the row and the las the column number.

Is it possible to produce the original raster file?


texts: Preventing labels from lower layers showing through layers that are on top - QGIS

texts:
I am creating a sequence of layers and layer groups - of geological units. Layers and groups higher up the sequence cover - partially or completely - lower layers or groups in the sequence, as they are made active, or visible. But I am finding that the labels for ALL active layers are showing, even through the actual visible layers - they are all polygon layers - which cover those previous layers. How can I get the labels to only show if that layer is on top? It is okay if only a part of a layer shows - as long as the label stays on that part of the layer that is visible.
I am using QGIS 3.16.1 Hannover on Windows 7.


texts:
I was wondering if anyone knows of a free WMS of imagery that I could use within QGIS.  The imagery does not have to be high resolution (Blue Marble even).  I would like to have a base imagery feed from a WMS.

Any help would be appreciated.

Thanks...


texts: Dynamically generating data for map tiles from Python function

texts:
I'm brand-new to GIS and relatively new to web services. I've found a ton of great resources on the web, but I'm having trouble addressing one specific issue:

One of the data layers on my planned map is dynamically generated in Python from some scattered and dynamic data (it is the "results" output of a computational weather-related model). I can generate this layer in chunks of arbitrary size and location, and output in any standard format (a numpy array, a .png, a geotiff (using GDAL), etc.), and I can expose a function such as:

def create_tile(width, height, ulc_lat, ulc_lon, lrc_lat, lrc_lon):

to hide the implementation details.

But, how do I get this layer to display in the map? The options that I can see right now are:

Serve as a WMS by creating a new Layer Type in a Python-based WMS server such as TileCache or TileLite ... this is probably what I'll try unless someone has a better idea.
Implement as a PythonDataSource in Mapnik. Problem: It's a new feature, requires custom compilation, and causes Mapnik to segfault even on the provided examples.
Somehow map filenames on disk (my_tile_z0_x3_y4) so that their contents is generated on-demand by the create_tile function above... I have no idea how this would be done, though.
Create the layer on disk ahead of time, and re-generate it often to keep the appearance of being "dynamic". (lots of disk space, needless CPU use, etc)

Does anybody have any suggestions? It seems like this must be a common problem, but I can't find any standard solution...

Thanks!
Rob


texts: Delaying load of ArcMap Addin extension?

texts:
For configuration reasons, I want to delay the load of my addin extension until ArcMap 10.2 has finished.

If there a property, flag, listener etc I can get a hook to do a check if arcmap is done?


texts:
I checked yesterday and today (2nd of December 2015).  The cartodb 'torque heat' disappeared from the wizard.

Is it something known or temporary?

Here is the information of what is (or should be) the 'torque heat' system:
http://blog.cartodb.com/introducing-heatmaps/


texts: I have just used the shape tool to create an elliptical buffer around a point. However after doing this, within my layers table I have the symbol seen within the figure below (seen in red box). Note, non of my other layers have any symbol next to them. Does anyone know what the symbol is and what does it mean? I ask this because when I use the Identify feature symbol it will only highlight the new elliptical buffer I just made and wont let me click on anything else.


texts: Fast solution to load and query shape files

texts:
I work on a .net core web API application collecting fishing events with latitudes and longitudes.

For each of these positions, I need to compute the corresponding Exclusive Economic Zone (EEZ), it's the national water surrounding each countries.

The official definition of these EEZ is available as shapefile on marineregions. There are few hundreds polygons.

For the moment I have stored the definition of these EEZ as polygon using the geography type in my SQL Server 2018 database, and I query with SQL the database for each position I collect.

I collect around 20,000 positions per day, and the performance is not great (takes 14 seconds for 500 positions).

So I'm looking for another solution.

What would be the fastest way to retrieve which polygon a point belongs to? Assuming I have shapefiles, but I can convert them to anything necessary. Also I need this solution to be free and if possible to offer an API so I can call it from my server.


texts: Creating/exporting multiple maps in one layout using ArcMap?

texts:
I have an air pollution data that consists monthly mean values of 12 monitoring stations. I made interpolation and want to export the 12 different map in one layout with one legend. Since I am a newbie with ArcMap, I wonder if it is possible.

I am using ArcGIS 10.1 for Desktop.

Edit: I mean - how can I automate apply symbology and export layout (for E.g. to JPEG) processes, so I wouldn't have to repeat same processes multiple times.


texts: What are the implications of modifying the Longitude of Origin?

texts:
I do a lot of GIS work with Alaska data, much of which comes from the State Geo-data Clearinghouse in EPSG: 3338 (NAD83 Alaska Albers), which has a longitude of false origin of -154.  Because Alaska is such a large area, making maps away from that longitude results in a skewed appearance where N, S, E, W lines are angled.  Here's an example from the northwest part of the state, using the default longitude value (-154) where the black lat/long lines are 68, -163:

In order to make such maps more pleasing to my audience, I often modify the CRS' longitude of false origin to match the longitude of the map center.  This can easily be done in either the QGIS Project Properties or ArcMap Data Frame Properties.  In the following example, I've changed the longitude of origin from the default -154 to match the map's -163:

I store the data, and do all analysis with the default CRS.  Modifying the latitude of origin is only done for final map production and display.

Obviously, the "non-square" effect occurs in other locations and with other CRSs.  I use Alaska for illustrative purposes.

AFAIK, there is no drawback to adjusting the longitude of origin for display.  To me, what I'm doing is akin to rotating a globe while I look at it, holding my eyes steady.  But still, is there something here that is cartographically improper?


texts:
I would like to calculate line sinuosity for line segments.

The sinuosity is defined as a ratio of the distance between start and end point along a line (line length) and the euclidean distance between both points.
A large sinuosity signifies a large detour between start and end point and is shown as the red line in the next image:

Caculating this measure is straightforward in GRASS GIS with:

v.to.db map=lines option=sinuous columns=sinuosity

However, v.to.db is not part of QGIS processing and cannot be called from within the QGIS Python shell.

Is there another option to calculate line sinuosity in PyQGIS?

Here some links to related questions:

Measuring Sinuosity of Walking Paths using QGIS?
Measuring straightness of a curve segment (represented as a polyline)


texts: How to directly get a Google API map based on a cell tower object

texts:
How do I pass a cell tower object to the Google Maps JavaScript API or Embed API to get a map for wherever the tower happens to be?

I realize that I can accomplish the task in two steps by first using the Geolocation API and then passing the location to the Maps API, but the extra round trip seems unnecessary.


texts: Can't select layer in the Road Graph plugin QGIS

texts:
I'm trying to use the Road Graph plugin in QGIS 2.2.0.

I have tried a few line shape files, but it won't let me select the layer (the option is greyed out). I've tried adding a 'SPEED' field and a 'Type' field but its still the same. The CRS is in the project format.

I can't see any other reason it wouldn't work?

Thanks,

Joe


texts: Split Lines at Intersection of Other Lines

texts:
I have two line layers (roads and contours in this example) and I would like to split the roads layer every time it crosses/intersects a feature on the contour layer. I can't seem to find a tool that does this, as the intersect tool is expecting a polygon for at least one of the layers. It also appears that GRASS' v.clean only works on a single layer. It seems like a simple task, but I'm at a loss!

EDIT (per comments below): Trying to get this to work in PostGIS now, and having some very strange results.  If I do a subset of my dataset (one road, one contour that crosses it), and use the following SQL, I get the expected result, with the road split where the contour crosses it:

insert into "road_split"
select
(ST_Dump(st_split(r.geom,c.geom))).geom
from "road" r, "contour" c

However, if I do the exact same thing with the complete dataset, no roads are split.  What am I doing wrong here?

I should note that all features in both the road and contour table are LINESTRINGS.


texts:
I'm trying to reproject raster layers in QGIS. At the moment they are in ESPG:4326 and I want to convert them to ESPG:54009. (Using Save As... and changing the CRS). However, this produces a layer of nan values.

Is there something obvious I am failing to do, to make it behave in this way?


texts: Action ran severall times when clicking on overlapping features

texts:
I have a layer of points and I associated a Python action to that layer which basically opens a url in the embedded QWebView.

This works well.

The issue I am facing is that if I zoom out and click on a region where several features are overlapping, the QWebView opens and reopens as many times as there are features stacked there.

Is there a way to prevent this? like:

Stop previous action when a new one is launched

or

Prevent the action to be ran for features that are not the first one visible


texts: Problem with weird values in DTM and not showing correct elevation

texts:
I am still new to Gis and don't know if I understand what is wrong.

I am working in a ocean survey program with point clouds. As seen on the picture the minimum depth is 1,92 and maximum 12,02. I first export the point cloud to a .dtm file. After that I export it as an tif (can't change anything here, only choose the option for which file format to save it in like tif).

The problem is when I look at it in Qgis there are 4 bands created and the maximum and minimum value are all weird (maybe I who can't interpret it). Shouldn't the value of elevation be stored in the cells as 1.92 min and 12.02 and not 248 and 0.


texts: Looking for data sources for real road speeds

texts:
A data question really, though there may be a useful geo-processing solution.

Context: I am  calculating "drivetime" catchment areas  using PgRouting and OSM / OSM2PO.  If my employer/client is to use this as information for decisions, then I think the base data needs to be good.

Problem: The fact is that in a large city say, London speed limits exceed traffic speeds due to congestion. Therefore, my polygons are not really valid (except at times of low activity).

Question: Thus, I'd greatly value any help (all ideas good) with how to update my OSM speeds on an area by area or road by road basis.  (E.g: Is there crowd-sourced data (from mobiles)?  Can I get access to traffic flow data? Has somebody else tackled this?).


texts:
I am attempting to calculate the distance between two latitude/longitude points.  I have a piece of code that mostly works that I yanked from this post but I do not really understand how it works.

Here is the code:

<?php
// POINT 1
$thisLat = deg2rad(44.638);
$thisLong = deg2rad(-63.587);

// POINT 2
$otherLat = deg2rad(44.644);
$otherLong = deg2rad(-63.911);

$MeanRadius = 6378 - 21 * sin($lat1);

$xa = (Cos($thisLat)) * (Cos($thisLong));
$ya = (Cos($thisLat)) * (Sin($thisLong));
$za = (Sin($thisLat));

$xb = (Cos($otherLat)) * (Cos($otherLong));
$yb = (Cos($otherLat)) * (Sin($otherLong));
$zb = (Sin($otherLat));

$distance = $MeanRadius * Acos($xa * $xb + $ya * $yb + $za * $zb);

echo $distance;
?>

I have a couple questions:

what are xa, ya, za?  I understand that they are points on a 3D cartesian plane but where are they relative to?  The center of the earth?
How does this cos($xa * $xb + $ya * $yb + $za * $zb) calculate the distance between the points?  I know that in 2D I would do this:

Pythagorean Theorem
distance^2 = b^2 + a^2
distance = sqr((y2-y1)^2 + (x2 - x1)^2)

How accurate will this be?  There was some discussion about that on the other page.  But I specifically want to use the distance to tell if users are within the something like 10m, 20m or 50m of each other.  Will I be able to do this with good accuracy?
What should I use for $MeanRadius?  Is that a reasonable value?  I think that that value assumes that the earth is a elipse.


texts: Combining Query Builder with Atlas Generation

texts:
I'm using QGIS 2.16.1 - Nodebo.

I have a layer containing 200 polygons which has been styled on a categorical variable. All these polygons - constituencies, belong to regions. There are 6 regions. I want to automate making region level maps.
[![enter image description here][1]][1]

Right now, I go to the query builder, and filter on a region and get the map like the one below and save it
[![enter image description here][2]][2]

How can I automate tell QGIS to go through all regions and filter on that region and save that region map? (automation is what I'm looking for - through a Python script or Atlas Generation)