texts: Two machines running the same arcpy based software - only one is leaking memory

texts:
Ok, so I'm running a rather large Python application based on arcpy on two separate machines - one of them is leaking memory and one is not. Both are running the same operation in the software.

The leaking machines is:

A nearly new Lenovo T420s
Intel Core i5-2540M processor
8 GB DDR 3 RAM
Windows 7 x64 Professional (non-domain joined)
128 GB Toshiba SSD
running ArcGIS 10.0 SP4 with pyodbc installed into the Python directory(only third party module installed)

The nonleaking machine is:

A 1.5 year old Dell Optiplex 980
Intel Core i7 870
8 GB DDR 3 RAM
Windows 7 x64 Professional (domain joined)
500 GB Western Digital 7200RPM HDD
running a brand new (today due to another issue with ArcGIS) intallation of ArcGIS 10.0 SP4 with pyodbc installed (also the only 3rd party module

I also have this same software running on other machines without leaking, but I haven't run this same operation recently or on SP4, so I'm leaving them out.

The software I'm running is here: https://bitbucket.org/nickrsan/pisces - The segment of it that's running is the mapping module (PISCES.mapping). Basically, it runs off the command line and is taking is input data from a Personal Geodatabase, generating layers based upon returned keys from stored queries, and adding those layers into template mxds and outputting them as PDFs/MXDs. All of the leaking portion uses arcpy.mapping and I suspect that's the culprit.

The nonleaking machine runs the software utilizing approximate 200 MB of RAM the whole time. The leaking machine starts around 150 MB, and slowly winds its way up to around 1.25 GB where arcpy calls (load layer, output PDF, etc) start failing.

Has anyone had any similar experiences? I know arcpy leaks memory sometimes and have had similar (but less severe) situations before. Any suggestions on how to track down a leak like this? Any best practices on the types of things I can do/avoid in arcpy to not have memory leaks? I've been using the gc module, but I think I'm not reading its output correctly. I'm relatively certain that it's not my code, but who knows.

Thanks in advance for your help.


texts: Get basic image/utf grid tile URLs from CartoDB visualization?

texts:
I'd like to add the layers from a CartoDB visualization to a Leaflet map directly (NOT through CartoDB.js).  My visualization has one image tile layer generated from data, and one utf grid tile layer generated from that same data, for interactivity.  How can I get from the big vizjson blob CartoDB produces to plain URL patterns like:

http://blah.cartocdn.com/blah/api/v1/map/{z}/{x}/{y}.png
http://blah.cartocdn.com/blah/api/v1/map/{z}/{x}/{y}.grid.json

This doesn't seem like it should be that hard, but the URLs that get used seem to be generated on the fly in some convoluted way.

The CartoDB.js API does have a getTiles method that should return those URLs, but it doesn't work.  It expects the layer options to include a "sublayers" property, which doesn't exist for me.  I gather those only exist for "anonymous" maps, rather than ones created using the CartoDB editor while logged in.

Any insight?


texts: Trying to publish a shapefile on GeoServer but the OpenLayers link on Layer preview can't find layer

texts:
I am having difficulties opening the OpenLayers link on the layer preview on GeoServer.

First I have trouble setting up a workspace on GeoServer. What is the workspace uri? I tried using one from the tutorial and from the demos present upon the installation of the server.

I then followed the various steps on the online documentation but in the end I get the download of this ugly XML file telling me the layer was not found.

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE ServiceExceptionReport SYSTEM "http://localhost:8080/geoserver /schemas/wms/1.1.1/WMS_exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1" >
<ServiceException code="LayerNotDefined" locator="layers">
      Could not find layer tutoWorkspace:final polygons
</ServiceException>
</ServiceExceptionReport>

Where am I possibly going wrong?

The namespace URI (which I still do not understand the criteria for its choice) or some other basic step.

I will be running the server offline.


texts: How to find boundary data for cities in Utah

texts: I'm working on a brand new project where I need boundary data for various cities in Utah. Salt Lake City was easy since their data is available for public download. The other cities? Not so much. Apparently this information is available as US Census data but I'm having a hard time parsing out the boundary data within. I'm more of a web developer so, I've never been in the position of tracking down the data itself. Annoyingly, a lot of cities in Utah have this information available as pdf download but I need the lat long city boundaries to build polygons in a web map. Can any of you more experienced cartographers help me?


texts: Substituting inline variables in ArcGIS ModelBuilder?

texts:
I want to iterate through a field of a feature class within a geodatabase in ArcGIS ModelBuilder, extracting each group to their own separate feature class within that same gdb.

I'd like to append the names of these groups to the end of the feature class name. As there are spaces and other invalid characters in many of the group names, I tried the Calculate Value tool to use Arc's ValidateTableName function to automatically replace the offending characters with underscores.

The model kind of works but I have an interesting anomaly where (as you can see above) an extra value of '1' has been added to the output 'Name' which serves as a precondition to the Calculate Value function. This results in '1' being added as an extra inline variable to the output feature names, causing a mismatch between the file names and the group upon which they were originally based.

How can I stop '1' from being added as an inline variable here?

Edit: interestingly if I remove the Calculate Values section the outputs are named correctly until it hits a group with an invalid character for geodatabases, at which point it falls over:


texts:
I work with LS8-data in a mountainous area where the illumination is poor especially in the winter months. My thought is to use either a shadow mask or apply topographic correction.

I will use all the available data for my research area over five years - aka apply a multi-temporal approach year-wise for a NDVI-time series.
Is there any way to implement e.g. the C-Correction in ENVI or what would be the best way to handle the topographic shadows?


texts: Converting Long/Lat into X/Y in the same coordinate system

texts:
I have a text file with two fields including Longitude and Latitude in degrees. I need to convert this into two fields that contain these values in X and Y. The source and destination coordinate system have to be the same: ESPG: 4236 - WGS 84.

I have tried the AttributeReprojector but it does not seem to work and the new X Attribute and Y attributes seem to be the same values as the original Longitude and Latitude.

Does anyone have an idea of how this could be resolved in FME?

I am showing a sample of the data below:

What I basically need is to get the X and Y corresponding to the Long/Lat in the projection WGS 84.


texts: Distorted Google Maps layer in pdf export from QGIS 2.0.1

texts: Export of a map from QGIS 2.0.1 Composer to pdf results in distortion (changed scale) of Google map layer. Why?


texts: Using OpenStreetMap Transport layer as xyz in OpenLayers

texts:
How can I use Openstreet Transport layer as baselayer in Openlayers?

var transport_layer =  new ol.layer.Tile({
    source: new ol.source.XYZ({
        attributions: '',
      url : 'https://a.tile.openstreetmap.org/${z}/${x}/${y}',

    })
  })


texts: QGIS does not seem to make use of spatial index when using GeoPackages

texts:
I have a QGIS project with several Geopackage-Layers.
The project is quite slow in comparison to a project with shapefile or database layers.
I started QGIS together with DebugView to find out what happens and recognised that every time I zoom in/out or select features QGIS makes a SQL select which queries every feature of the layers (Select * from...):

[248] ..\..\..\..\..\..\src\providers\ogr\qgsogrprovider.cpp(3504) : (QgsOgrProviderUtils::setSubsetString) [7000ms] SQL: SELECT *, fid as orig_ogc_fid FROM "lines" WHERE layer IN ('Verm. Randsteine und Belag') AND space=0 AND block=-1

I opened the file with SQLiteDatabaseBrowser to check if a spatial index exists and found a table "rtree_lines_geom" which seems to be used as a index:

CREATE VIRTUAL TABLE "rtree_lines_geom" USING rtree(id, minx, maxx, miny, maxy)

Is this a bug that QGIS always selects all features from the table instead of using an index?


texts: What does "The ArcGIS Server Object Manager service hung on starting" indicate?

texts:
I often get this error message in the Windows 2008 R2 server event logs:
The ArcGIS Server Object Manager service hung on starting. (EventID: 7022)

Have tried searching the cause for this error but noluck yet.
Any ideas?

S/W specs:

ArcGIS Server Advaced Enterprise 10 SP4

Regards
Ujjwal


texts: Display XY using ModelBuilder in ArcGIS Desktop?

texts:
My final results from a model builder is a .dbf file with X and Y coordinates as shown in the image.

Is there any way to plot these coordinates on ArcMap using ModelBuilder?


texts: Adding feature to GeoServer vector layer using WFS-T in absence of fid index?

texts:
I want to add a feature to GeoServer vector layer with WFS-T.

When I use saveStrategy.save(), an error returned that said:

  java.io.IOException: current fid index is null, next must be called before write().

How can I add a geometric feature to GeoServer vector layer without fill feature fid?


texts: In PostGIS how can I do an ST_Intersection() search with a GEOMETRYCOLLECTION Geography field?

texts:
I have a Postgres table "files" with a lat/lon location field of type geography(GeometryCollection,4326). Right now the field only contains lat/lon polygons. I have succeeded in creating a gist index for the field and in doing a bounding box search with a search polygon:

select fileid, file_name, ST_AsText(location)
from files
where location && ST_MakeEnvelope (-81.0, 22.0, -76.0, 27.0);

But, I need to do a fully intersection search but I hit the restriction that ST_Intersection() does not accept a GeometryCollection. Maybe I could use ST_Distance() since I only need to know if a search area intersects a location area, not what the intersection area is. But I don't see any way to specify a spherical earth calculation, and that is probably all I need most of the time. I have been thinking about changing the field to a multipolygon since that would probably work.

It looks like ST_Distance() will work for what I need to do. It accepts a GeometryCollection and actually will do spherical earth calculations.


texts: Generate map coordinates to be on streets

texts:
I have a PHP program that generates a certain number of random coordinates on a map in a specified radius. I want check which points are positioned in a user-accessible area such as a street; not on buildings.

I was thinking of generating a query for directions, getting the highlighted map and making sure the points are positioned on the lines.

Any ideas?


texts: Using the watershed tool with two pour points

texts:
I have two pour points in a stream network created using the flow accumulation methodology in ArcMap.

When I use the watershed tool with this pour point shapefile, I get a watershed between these two points. What exactly is this watershed?

I've tried searching on the ESRI site, but I can't seem to pinpoint the exact answer.

Edit: I don't want two catchments here. I want an area of everything that flows from the "upper pour point" to the "lower point." Picture included.


texts:
Working on a Python script that I will run inside of ArcGIS 10.5 that will export PNG maps. Each map will be a different state and for every state there will be 5 different physician specialty types, resulting in about 250 total maps. My current version of the script can export one type of physician for all states, but when I move on to the next physician type, I have to manually enter the next one into my script and run again.

I can loop through all 50 states, but I am not finding success in looping through all 5 physician types in a single go. I suspect it's part Python and part data set-up, but I was wondering if anyone had pointers on how to effectively loop within a loop.

So right now, it loops through a dictionary list of the states. How do I include a sub-loop list and be able to use it in a for loop in the script?

Code is below; all of the 'spec%' items are the physician types and in my MXD, I have separate layers for each point layer.

for state, fname, spec1, spec2, spec3, spec4, spec5, spec6, spec7, spec8, spec9, spec10, spec11, spec12, spec13, spec14, spec15, spec16, spec17, spec18, spec19, spec20, spec21, spec22, spec23, spec24, spec25, spec26, spec27, spec28, spec29 in st_list_tup:

lyr3.visible  = False
lyr4.visible  = False
lyr6.visible  = False
lyr7.visible  = False
lyr8.visible  = False
lyr9.visible  = False
lyr10.visible = False
lyr11.visible = False
lyr12.visible = False
lyr13.visible = False

# CREATE SELECTION LAYER AND WHERE CLAUSE
arcpy.SelectLayerByAttribute_management(layer,"NEW_SELECTION", st_abbr.format(state))

# ZOOM TO STATE SELECTION THEN CLEAR SELECTION
df.zoomToSelectedFeatures()

# APPLY SYMBOLOGY
arcpy.ApplySymbologyFromLayer_management (lyr6, symlyr3)

# CLEAR SELCTION
arcpy.SelectLayerByAttribute_management(layer, "CLEAR_SELECTION")

# DEFINITION QUERIES
lyr.definitionQuery = st_abbr.format(state)
lyr5.definitionQuery = count_lyr.format(fname)

for layer in arcpy.mapping.ListLayers(mxd):
    if layer.name == "Facial Surgeons":
        layer.visible = True
        layer.definitionQuery = FSurg_lyr.format(state)

print lyr.definitionQuery, lyr3.definitionQuery, lyr4.definitionQuery

# UPDATE MAP ELEMENTS
for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
    if elm.name == "Title":
        elm.text = "{}".format(spec16)

for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
    if elm.name == "Subtitle":
        elm.text = "{}".format(fname)


texts:
I am Trying to recreate this working python function in Javascript:
task = ee.batch.Export.image.toCloudStorage(
    image = image,
    description = 'test',
    bucket = BUCKET,
    fileNamePrefix = FOLDER + '/' + 'test',
    region = aoi.getInfo()['coordinates'],
    crs = PROJ["crs"],
    crsTransform = PROJ["transform"],
    fileFormat = 'TFRecord',
    maxPixels = 5e10,
    formatOptions = {
      'patchDimensions': KERNEL_SHAPE,
      'kernelSize': kernel_buffer,
      'compressed': True,
      'maxFileSize': 104857600
    }
    )
task.start()

The goal is to send the image to google drive instead of the storage bucket
When running
Map.addLayer(image,[],"image")
print("image.projection().getInfo()",image.projection().getInfo())
var proj = image.projection().getInfo()
var KERNEL_SIZE = 1
var kernel_buffer = [ee.Number((KERNEL_SIZE+1)/2).int(), ee.Number((KERNEL_SIZE+1)/2).int()]
print("kernel_buffer",kernel_buffer)
var aoi_name = "test"
var mode = "test"

Export.image.toDrive({
        image : image,
        description : 'Whole_Tfrecords_'+aoi_name,
        folder: "GEEExports",
    fileNamePrefix:"Europe_"+aoi_name+"WHOLE_TF_record_"+mode,
        region : aoi,
        crs : proj.crs,
    crsTransform : proj.transform,
        fileFormat : 'TFRecord',
        maxPixels : 5e10,
        formatOptions : {
          'patchDimensions': kernel_buffer,
          'kernelSize': kernel_buffer,
          'compressed': true,
          'maxFileSize': 104857600
        }
  }
        )

I get this Error.
Unable to construct grid from dimensions: NaN,NaN,NaN,NaN,NaN,NaN

The Documentation states:
region (Geometry.LinearRing|Geometry.Polygon|String, optional):
A LinearRing, Polygon, or coordinates representing region to export. These may be specified as the Geometry objects or coordinates serialized as a string. If not specified, the region defaults to the viewport at the time of invocation.
Playing with the region and using the serialisation function did not yield results.
neither works
var = ee.List(ee.Geometry(aoi).coordinates()).serialize()
var = ee.Geometry(aoi)
var = ee.List(ee.Geometry(aoi).coordinates())
var = ee.Geometry(aoi).coordinates()
var = aoi.getInfo().coordinates.serialize()
var = aoi.getInfo().coordinates
var = ee.Geometry(aoi).serialize()
var = ee.List(ee.Geometry(aoi).getInfo().coordinates).serialize()
Link to GEE Script


texts: what is the right way to delete all features from featureLayer

texts: I have a featureLayer and I don't want to kill it and create another one. I want to delete all features from the current featureLayer. I found that I can use applyEdits method but in order to use it I need to have all features I want to delete. Whilst I don't want to store them separately. Do you know a better way to empty my featureLayer when using the ArcGIS Javascript API?


texts: Why am I unable to create overviews in QGIS despite having disk space

texts:
I have 120 GeoTIFF files, totalling around 250GB of data.

Files names are from 1.tif to 120.tif.

I tried to create overviews for all the files in one go usinggdaladdo in terminal. While overviews got created till 50.tif, for the rest of the files I got the no more space on disk error. I checked and indeed my drive partition (264GB) was full.

I moved the rest of the files to another partition on my disk, which had over 200 GB available.

Now when I create overviews using gdaladdo for these relocated files in their new location, I get the same error. This time I don't create the overviews using gdaladdo in terminal. Instead, I use QGIS. (I know QGIS also uses gdaladdo only internally, but I thought I should mention this.)

ERROR 1: /media/nancy/G/51.tif:Failed to allocate memory for to read TIFF directory (0 elements of 12 bytes each)

ERROR 1: TIFFReadDirectory:Failed to read directory at offset 1195893042

0...10...20...30...40...50...60...70...80...90...100 - done.

Execution completed in 2.04 seconds
Results:
{'OUTPUT': '/media/nancy/G/51.tif'}

Also, on loading 51.tif in QGIS and right clicking, I get a cross on all the resolutions of the overviews.

This is not the case for some other rasters that I have in my /home directory. This same window for them shows all the overview resolutions possible to be made.

What is wrong?


texts:
May i know whether QGIS supports the following map formats?
(1) S-57 Electronic Navigation Charts (Naval Maps)
(2) Army/IAF Maps in DGN
(3) DTEDO0, DTEDO1 and DTEDO2
(4) GeoTIFF
(5) DTED


texts:
In a weeks time I should be able to transfer the GIS data provided in the XML format to PostgreSQL database. I have prepared the PostgreSQL database (empty) and looking at Matlab scripts. (I am on ubuntu platform).

Matlab is our existing platform.The work is passed on to me by another person who used Mapping tools on MATLAB to export the data to viewers like ArcGIS/Google Earth.This was created for a prototype. (However, in the final product I will be exporting PostgreSQL to QGIS to view the data).

I am pretty new to the GIS data/table.Could you suggest me the steps in-detail to have the data in PostgreSQL?


texts: Searching multiple geodatabases for feature classes and merging the results

texts:
I've got some Sea Level Rise data from NOAA (https://coast.noaa.gov/slrdata/) stored in 66 different geodatabases. Each geodatabase represents a state or district and contains a set of feature classes of sea level rise at different depth (0 feet, 1 feet,...). As is shown in these pictures:

My aim is to merge all the features that have the same water depth, for example, 0 feet, to a single feature class to create a sea level rise map of the US at different water depths. So I used Arcpy to find polygon feature classes that end with "slr_0ft" through these geodatabases.  And then merge the results to one single feature class. Here's the code I used:
import arcpy, os
workspace = "G:\DirectStudy\SLR\gdb" #My workspace that contains 66 GDBs
search = "slr_0ft"
feature_classes = []
for dirpath, dirnames, filenames in arcpy.da.Walk(workspace,
                                                  datatype="FeatureClass",
                                                  type="ANY"):
    for fname in filenames:
        # search for string in string to eliminate the need for exact filenames
        if search.upper() in fname.upper():
            feature_classes.append(os.path.join(dirpath, fname))
    for fc in feature_classes:
         print fc

arcpy.Merge_management(feature_classes, "G:\DirectStudy\SLR\Merge_0ft")
print "Done!!"

The code was successful in printing paths and filenames (see below). But the output file from the merge function is an empty shapefile. I checked it with ArcMap and found it contains no polygon. Also, the script didn't give me any kind of error or warning. I guess the mistake should be in the arcpy.Merge_management() function. How can I fix my script?

I'm not a native speaker so I tried best to explain my question clearly.


texts:
I know about using ogr2ogr, and also about the nodejs module osm_geojson. However, I would like to convert GeoJSON to OSM XML from Python. Is there a library which has this capability built in? Or would I have to write custom code based on the OSM XML specification?
I am processing the GeoJSON using geopandas, so exporting directly from geopandas to OSM would also work.

I have looked at the GDAL/OGR Python library.
It may be able to do what I want, but the documentation is impenetrable. The Cookbook has something similar to what I am looking for, but I can't quite seem to piece it together.

My code currently looks like this:

from osgeo import ogr
import os

# Read GeoJSON file
geojson_file = 'file.geojson'
GeoJSON_driver = ogr.GetDriverByName('GeoJSON')
data = GeoJSON_driver.Open(geojson_file, 0)
layer = data.GetLayer()

# OSM file to write
osm_file = 'file.osm'
OSM_driver = ogr.GetDriverByName('OSM')

# Remove output file if it already exists
if os.path.exists(osm_file):
    OSM_driver.DeleteDataSource(osm_file)

# Export data to OSM file
# This is where I'm stuck


texts: can wms from geoserver serve as utfgrid data?

texts:
i'm confused on how to make interactive map displayed. i used geoserver from opengeo suite and then i use openlayers to display the map on my site. I have been able to display the attribute by using getfeatureinfo function. The problem is that i need to be able to use utfgrid like this example http://openlayers.org/dev/examples/utfgrid.html to show one of my attribute("location").
How can i do it? i've been searching all day about this, all information leads to use tilemill/stache etc rather than geoserver.
any help would be appreciated. thanks


texts: Concatenating coordinates inside geometryFunction

texts:
While looking at this example, I came upon the following code (taken verbatim by using 'view page source'):
interaction: new ol.interaction.Draw({
      type: 'Polygon',
      source: vector.getSource(),
      // Count inserted points
      geometryFunction: function(coordinates, geometry) {
        this.nbpts = coordinates[0].length;
        if (geometry) geometry.setCoordinates([coordinates[0].concat([coordinates[0][0]])]);
        else geometry = new ol.geom.Polygon(coordinates);
        return geometry;
      }
    }),

I understand what the line starting with this.nbpts =  does in that particular example.
What I am most puzzled is the following line:
if (geometry) geometry.setCoordinates([coordinates[0].concat([coordinates[0][0]])]);

… in particular the weird concatenation it uses at the end. I experimented with replacing the above line with the much simpler:
if (geometry) geometry.setCoordinates(coordinates);

… and it works equally well, so I don't understand what the concat achieves.


texts:
I have the following code which batch saves rasters in qgis using the python terminal:

myDir = 'd:/output/'
layers = iface.legendInterface().layers()
pipe = QgsRasterPipe()
for layer in layers:
    extent = layer.extent()
    width, height = layer.width(), layer.height()
    renderer = layer.renderer()
    provider=layer.dataProvider()
    crs = layer.crs().toWkt()
    pipe.set(provider.clone())
    pipe.set(renderer.clone())
    file_writer = QgsRasterFileWriter(myDir + layer.name() + ".tif")
    file_writer.writeRaster(pipe,
                            width,
                            height,
                            extent,
                            layer.crs())

The tifs that are saved are uncompressed, does anyone know how to turn on compression like LZW?


texts:
I have a map with a lot of labels with the value 0. I want to hide or remove them somehow and came up with the following rulebased label:

if( floor("area") = 0, " ",floor("area")||'m²')

Unfortunately, now every label is hidden. But I don't see the error, it should work shouldn't it.

Perhaps the idea with the if question is not good.


texts: Black image after clipping with ArcGIS for Desktop?

texts:
I am trying to clip my tif image but i got it readable only in arcgis, and when i open the folder where it is recorded i found it black and cannot open it or read it with Autocad ..
what do you think about ?

I am using ArcGIS 10.1 and made the commande Spatial analys ==> Extraction ==> Extract by Mask

I tried also  Data Management tools ==> Raster ==> Raster processing ==> Clip

But the both give the same result.


texts:
I have a scanned image (tiff) of a historic map, with no information of the used coordinate system, which I want to georeference, to finally be part of a WebMap application with Leaflet.

As far as I know, for WebMaps the coordinate reference system is the "WGS84 Web Mercator (Auxiliary Sphere)". Due to that, I thought it would be best to georeference it directly to the used coordinate system.

My approach would be to georeference the historic image using an ArcGIS basemap.

My constraints are that ArcGIS basemaps do not have the best resolution, and I am worrying about the accuracy of georeferencing. What I could do is to include some vectors of topographic features (rivers, streets, lakes) to increase the accuracy, but I am not sure if this is a valid way to do.

Is it recommended to georeference based on ArcGIS basemaps or could you suggest a more appropriate approach?

The historic map is a topographic map, it shows landscape features and in addition old paths. These paths should be digitized, after the map is georeferenced, to create a WebMap with Leaflet showing users where to walk on "old paths".

I tried to georeference by mountain ranges, that could be seen on the historic map as well as on the basemap (although I am not sure in which quality these are represented on the basemap) and by river crossings which still exist. Just to clarify, I did not use the course of the river itself, but only intersecting points of two or multiple streams that are still valid.

I am georeferencing using ArcGIS Desktop.


texts:
I have a grid layer (lidar_cortina) with "NAME" and "percorso"(path) coloumns:

I'm trying to write an QGIS action that loads file .asc from a grid layer.

qgis.utils.iface.addRasterLayer("[%percorso%]", "[%NAME%]")

but It doesn't work and the message "the command can not be executed" appears

If I use the console loading one .asc it works well

qgis.utils.iface.addRasterLayer("//rvphnas02pw/Remote_Sensing/Lidar/2017_Lidar_Cortina/SECONDA_CONSEGNA/Lidar/DTM/P740451578_0101_100_DTM.asc", "P740451578_0101_100")

Any suggestion.

UPDATE

If I use this code, it works well:

qgis.utils.iface.addRasterLayer('//rvphnas02pw/Remote_Sensing/Lidar/2017_Lidar_Cortina/SECONDA_CONSEGNA/Lidar/DTM/[%NAME%]_DTM.asc', '[%NAME%]')

Using the code of: Developing QGIS Action which loads a raster
I have write an other qgis action to open .asc file from a grid layer with "NAME" column:

from PyQt4.QtCore import QFileInfo,QSettings
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem

vl = QgsMapLayerRegistry.instance().mapLayersByName('lidar_cortina')[0]
qgis.utils.iface.setActiveLayer(vl)
root = QgsProject.instance().layerTreeRoot()
group_name = "lidar"
group = root.findGroup(group_name)
if group == None:
    group = root.addGroup("lidar")
else:
    pass
fileName ='D:/GIS/Dati/lidar/[% NAME %]_DTM.asc'
fileInfo = QFileInfo(fileName)
baseName = '[% NAME %]'
qgis.utils.iface.addRasterLayer(fileName, baseName)
rlayer = QgsRasterLayer(fileName, baseName)
crs = QgsCoordinateReferenceSystem()
crs.createFromSrid(32632)
rlayer.setCrs(crs)
QgsMapLayerRegistry.instance().addMapLayer(rlayer, False)
group.insertChildNode(-1, QgsLayerTreeLayer(rlayer))

but It doesn't work.

Using python console to open one .asc file with the same code of error message, It works:

from PyQt4.QtCore import QFileInfo,QSettings
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem

vl = QgsMapLayerRegistry.instance().mapLayersByName('lidar_cortina_prova')[0]
qgis.utils.iface.setActiveLayer(vl)
root = QgsProject.instance().layerTreeRoot()
group_name = "lidar"
group = root.findGroup(group_name)
if group == None:
    group = root.addGroup("lidar")
else:
    pass

fileName ='D:/GIS/Dati/lidar/P740451578_0101_100_DTM.asc'
baseName='P740451578_0101_100'
fileInfo = QFileInfo(fileName)
rlayer = QgsRasterLayer(fileName, baseName)
crs = QgsCoordinateReferenceSystem()
crs.createFromSrid(32632)
rlayer.setCrs(crs)
QgsMapLayerRegistry.instance().addMapLayer(rlayer, False)
group.insertChildNode(-1, QgsLayerTreeLayer(rlayer))

I can't understand
Any suggestion is welcome.

UPDATE 2

changing the path with the suggestion in comment of @lrssvt, I obtain the same error:


texts:
I am trying to open a Jp2 file in QGIS 3.4.5 on a macOS 10.14.5 Mojave.

I am new to QGIS and have been looking through forums that discuss GRASS and GDAL but I don't understand the connection to QGIS.

DETAILS:
- When I "drag and drop" I get the dreaded beachball and QGIS remains unresponsive.
- The file is a statewide raster of magnetic data, roughly 120 MB. It also won't open in preview, with the following error "It may be damaged or use a file format that Preview doesn’t recognise".
- The GDAL version is 2.4 which I downloaded and installed from www.kyngchaos.com.

Can you give me a simple explanation?


texts: Assign elevation data to GPX from LIDAR data

texts:
I am totally new to this, so sorry for the basic question.

I have some GPX files which have elevation data, but I am aiming for greater accuracy. Currently I have used OS benchmarks to correct the barometric elevation data. Online digital elevation models (ASTER, SRTM, Strava, etc) do not provide the degree of precision or accuracy that I was seeking and perform more poorly than the barometric data.

I have downloaded QGIS, imported a shape file as a layer from the Environment Agency (UK) and imported a GPX route. Unfortunately the shape file appears, as I suspected, to be just the areas covered by the LIDAR and not the actual elevation data.

Does anyone have a dummies guide to help me please? It seems I am stuck at a most basic level.

Also, later I would actually like to be able to draw a GPX route and get elevation data from the LIDAR dataset and export as GPX. Is this possible do you think?

I would like to be able to access much of the data linked to on this page: https://www.researchgate.net/post/Which_EU_Countries_provide_free_access_to_LIDAR_data_files_in_2017 and more if available.

EDIT:
I have got to here:

But I am really struggling to see how to get the elevation assigned to the GPX despite the link in the answer below. (and to export gpx when that is done). Any help would be most appreciated please.


texts: Trouble with polygon to polygon "Join attributes by location"

texts: I have two polygon vector layers - one of census tracts and one of water bodies. I would like to join a particular attribute from the water bodies table to each census tract, with the output being the sum of this attribute for all water bodies that border or intersect that census tract. I'm using the "Join attributes by location" tool, with summary of intersecting features checked. First, it tells me that there is a CRS error, which there is not - both layers are in the same CRS (EPSG:3310 - NAD83 / California Albers). Then, although it does not give me any errors, the output file only has null values for all the water body layer's features. The polygons in the output file are correct, however. Is there something I'm doing wrong here? Using version 1.8. Thanks!


texts: Couldn't load plugin PRISM due an error when calling its classFactory() method in QGIS1.8 inWin8

texts:
When I load a plugin in QGIS ,met some problems ：

I had tried copy other Python27's files to QGIS's Python 27 file but it still does't work.

Couldn't load plugin PRISM due an error when calling its classFactory() method

Traceback (most recent call last):
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 164, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins\PRISM__init__.py", line 25, in classFactory
    from prism import PRISM
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins\PRISM\prism.py", line 15, in
    import numpy
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\numpy__init__.py", line 137, in
    import add_newdocs
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\numpy\add_newdocs.py", line 9, in
    from numpy.lib import add_newdoc
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\numpy\lib__init__.py", line 4, in
    from type_check import *
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\numpy\lib\type_check.py", line 8, in
    import numpy.core.numeric as _nx
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\numpy\core__init__.py", line 5, in
    import multiarray
  File "C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 309, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: DLL load failed: %1 ���O���Ī� Win32 ���ε{���C

Python 版本:
2.7.5 (default, May 15 2013, 22:43:36) [MSC v.1500 32 bit (Intel)]

QGIS 版本:
1.8.0-Lisboa Lisboa, 6416f38

Python 路徑: ['C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users/corbies/.qgis//python', 'C:/Users/corbies/.qgis//python/plugins', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\Program Files\Intergraph\ERDAS IMAGINE 2014\usr\lib\Win32Release\python', 'C:\WINDOWS\SYSTEM32\python27.zip', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\DLLs', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\plat-win', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\lib-tk', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\bin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages', 'C:\Program Files (x86)\ArcGIS\Desktop10.2\bin64', 'C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy', 'C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Scripts', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\python\plugins\fTools\tools']

I think that the plugin designed for QGIS 1.8-1.99


texts: In QGIS 3.4 is there a tool similar to the Extract Layer Extent but instead of a bounding box, it creates polygons? I have groups of lines representing a network of roads and I want to show the missing sections but creating polygons around the groups. How can I achieve this?


texts: Create points along 3D stream where elevation changes by 1ft

texts:
Using ArcPy (Python 3), I am trying to create points along a line where the elevation changes by 1 ft between vertices. I can successfully convert my Stream Polyline feature class to 3D and create a Point feature class where the elevation changes every 1 ft. I now need to add the Z elevation from the vertex to the Point feature class. Can I do by using an Update Cursor within the Insert cursor? I'm still pretty new with Python, what I remember from class is that it's bad practice to create nested cursors? What other approach could I use to do this?

Update: Was able to figure it out. If I create point feature class that is Z aware the point vertices will have the elevations. Then I can use the AddZInformation tool to transfer those values to the Point's attribute table.

Update: Now I am trying to create points every 1000 feet between my existing points if no points currently exist. How can I do this? I'm thinking the Generate Points Along Line tool might work. https://pro.arcgis.com/en/pro-app/tool-reference/data-management/generate-points-along-lines.htm

import arcpy
import os
from arcpy import env

arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("3D")
arcpy.CheckOutExtension("Spatial")
env.workspace = r"C:\GIS\Tools\BFE_Tool\BFE_Tool.gdb"

streams = 'FY16_Streams_Diss' # Streams must be dissolved by stream name
outGDB = r'C:\GIS\Tools\BFE_Tool\BFE_Tool.gdb'

print ('Interpolating shape...')
stream3D = arcpy.InterpolateShape_3d(r'C:\GIS\Tools\BFE_Tool\Test_Data\wsg100x', streams, 'FY16_Streams3D')
print ('Interpolate shape done!')

print ('Creating BFE Points...')
try:

    myLayer = stream3D
    prev = -1
    heightDiff = 1
    pntList = list()
    with arcpy.da.SearchCursor(myLayer,["SHAPE@"]) as cursor:
        for row in cursor:
            geom = row[0][0] # A whole polyline
            for i,point in enumerate(geom):
                if prev == -1:
                    prev = point.Z
                else:
                    if abs(point.Z - prev) > heightDiff:
                        pntList.append(geom[i-1])
                        pntList.append(geom[i])
                    prev = point.Z

    if len(pntList) > 0:
        print('Creating Empty point feature class...')
        outFC = arcpy.CreateFeatureclass_management(outGDB, 'BFE_Points', 'POINT','', '','ENABLED') # Create Empty Feature Class to store Points

        # print ('Adding Elevation field to Points...') # Add field in Points to store Water Surface Grid Elevations
        # arcpy.AddField_management(outFC, 'WS_Elev', 'FLOAT')

        with arcpy.da.InsertCursor(outFC, ["SHAPE@"]) as cursor:
            for point in pntList:
                print(point.Z)
                cursor.insertRow([point])
                print ('Point Created')

    print ('Adding Elevation to Attribute table...') # Add field in Points to store Water Surface Grid Elevations
    arcpy.AddZInformation_3d(outFC, 'Z')

except arcpy.ExecuteError:
    # Geoprocessor threw an error
    print(arcpy.GetMessages(2))

print ('Finished!')


texts: Shapefile and Database file too large for Google Earth/Accumap

texts:
I am completely new to GIS and I had hopefully a simple question regarding how to reduce the file size of my shapefile and associated database file.  The files sizes are:

Shapefile (.shp file format) is 50 mb
Database (.dbf file format) is 50 mb

The attribute table information I have is:

12 fields
length averaging 10-20 with precision averaging 0-1
In total there is 417960 lines in the attribute table for each one of the 12 fields

With my data I also have a .prj file and .shx file. I have been trying to use QGIS 2.18.13 with GRASS 7.2.1 to convert the files into KML and also to change into NTD83 for use in IHS Accumap.

When I try to import the files into Google Earth the error message states there is over 2500 layers and even when viewing in restricted view it slows and crashes Google Earth.  And IHS Accumap is very slow as well.


texts:
I'm trying to retrieve the following pixel statistics:

Proportion of zero value pixels:

[nº of zero value pixels / (total nº of pixels - nº of NA pixels)]
Proportion of pixels smaller than 3rd quartile (= 0.75):

[nº of pixels < 3rd quartile / (total nº of pixels - nº of NA pixels)]

Here is a reproducible example:

library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f) #r is the object
plot(r)
summary(r)


texts: How to include maps in Jupyter notebook HTML export?

texts:
I'm running the 1.4.2 version of the ArcGIS Python API. In a Jupyter Notebook, I create a map like:

lowes = "Arkansas USA Lowe's"
lowes_map = gis.map("AR, US", 6)
lowes_map

...and get my output. I then geocode some points and plot them:

All is great there, but when I export to HTML, my map isn't there (this particular notebook has maybe 8 maps and none of them get exported to HTML):

If I display Portal items in my notebook, they get included in the HTML export. Also, all the code samples on the Python API docs have maps in them (such as here), so I'm thinking I'm just missing something?

Anyone know how to get maps to render in HTML export?


texts: Does anyone know of any GIS tool that reorders polygons so that the path through all polygons is minimized, i.e. it solves the travelling salesman problem? Here's an example of just how poorly optimized a shapefile can be:


texts: pgRouting on RHEL6: issues with driving distance library

texts:
I installed pgrouting on RHEL6 via yum and the ELGIS 6 repo (the only place where I could find pgrouting for RHEL6).

The problem is that the librouting_dd.so library (for driving distance) is missing. In other words, I have librouting.so, librouting_tsp.so, but not librouting_dd.so. Does anyone know why this is or how to fix it? I've included the yum output below.

I also tried compiling pgrouting from source (with the -DWITH_DD=ON flag so librouting_dd.so would be built) but this resulted in a waterfall of dependency issues (e.g., it needs CGAL-devel, which in turn needs things like lapack-devel and blas-devel which aren't installing because of other dependency issues, etc., etc.).

Any advice would be greatly appreciated.

# yum install pgrouting
Loaded plugins: rhnplugin
This system is receiving updates from RHN Classic or RHN Satellite.
elgis                                                                                | 3.0 kB     00:00
elgis-plus                                                                                | 1.3 kB     00:00
elgis-source                                                                                | 3.0 kB     00:00
Setting up Install Process
Resolving Dependencies
--> Running transaction check
---> Package pgrouting.x86_64 0:1.05-1_0.el6 will be installed
--> Finished Dependency Resolution

Dependencies Resolved

============================================================================================
 Package         Arch      Version         Repository                      Size
============================================================================================
Installing:
 pgrouting      x86_64     1.05-1_0.el6    elgis                            66 k

Transaction Summary
============================================================================================
Install       1 Package(s)

Total download size: 66 k
Installed size: 226 k
Is this ok [y/N]: y
Downloading Packages:
pgrouting-1.05-1_0.el6.x86_64.rpm                                                                                |  66 kB     00:00
Running rpm_check_debug
Running Transaction Test
Transaction Test Succeeded
Running Transaction
  Installing : pgrouting-1.05-1_0.el6.x86_64                                                                                1/1
  Verifying  : pgrouting-1.05-1_0.el6.x86_64                                                                                1/1

Installed:
  pgrouting.x86_64 0:1.05-1_0.el6

Complete!
[root@vm1 yum.repos.d]# rpm -sql pgrouting
normal        /usr/lib64/pgsql/librouting.so
normal        /usr/lib64/pgsql/librouting_tsp.so
normal        /usr/share/doc/pgrouting-1.05
normal        /usr/share/doc/pgrouting-1.05/README.routing
normal        /usr/share/postlbs
normal        /usr/share/postlbs/matching.sql
normal        /usr/share/postlbs/routing_core.sql
normal        /usr/share/postlbs/routing_core_wrappers.sql
normal        /usr/share/postlbs/routing_topology.sql
normal        /usr/share/postlbs/routing_tsp.sql
normal        /usr/share/postlbs/routing_tsp_wrappers.sql


texts: Why is the map frame suddenly empty in QGIS-print composer?

texts:
Often when working in QGIS  print-composer (2.14 and prior), at a certain point the map content disappears and only the empty map frame remains. This happens mostly after closing and re-opening my project or randomly after I hit the "refresh preview" button.

My solution so far was to just add new mapframes (over and over again) but this is annoying.

What is the reason behind my map-preview disappearing and what to do about it? The screenshot shows the problem and the general adjustment I am using.


texts: Multiprocess Spatial Join with Spatial Index

texts:
I'm a newbie in implementing GIS on Python, but the lack of any GIS software at the office has forced me to learn about it.

I have a huge grid (a multipolygon) that I'm currently filling out with counts of points using spatial index. When using a subset of the grid the code works just fine but when I use anything bigger than that then it takes forever.

This is the code that I'm currently using

# build the r-tree index
sindex = points.sindex

# find the points that intersect with each subpolygon and add them to points_within_geometry
points_within_grid = pd.DataFrame()

# grid-level dataseet
grid_data = pd.DataFrame()

i=1
for poly in grid:
    clear_output(wait=True)

    # buffer by the <1 micron dist to account for any space lost in the quadrat cutting
    # otherwise may miss point(s) that lay directly on quadrat line
    poly = poly.buffer(1e-14).buffer(0)

    # find approximate matches with r-tree, then precise matches from those approximate ones
    possible_matches_index = list(sindex.intersection(poly.bounds))
    possible_matches = points.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.intersects(poly)]

    # count number of firms within polygon
    x=pd.DataFrame({'grid_index': i,
                    'geometry':[Polygon(poly.envelope)],
                    'n_firms': len(precise_matches)})

    for poi in range(1,9):
        x['sic_'+str(poi)]=len(precise_matches[precise_matches.sic1==poi])

    grid_data=grid_data.append(x)

    # points within the grid
    points_within_grid = points_within_grid.append(precise_matches)
    print("Current Progress: ", np.round(i/len(grid)*100,0),"%")

    i=i+1

# drop duplicate points, if buffered poly caused an overlap on point(s) that lay directly on a quadrat line
points_within_grid = points_within_grid.drop_duplicates(subset=['x', 'y','poi_id'])
points_outside_grid = points[~points.poi_id.isin(points_within_grid.poi_id)]

So, I've been thinking about using multiprocessing.Pool to make it faster. I've been reading online but don't really know how to do it. The first step would be to define a function that can be introduced to the Pool.

So far I have this function,

def grid_inter(poly)
    # buffer by the <1 micron dist to account for any space lost in the quadrat cutting
    # otherwise may miss point(s) that lay directly on quadrat line
    poly = poly.buffer(1e-14).buffer(0)

    # find approximate matches with r-tree, then precise matches from those approximate ones
    possible_matches_index = list(sindex.intersection(poly.bounds))
    possible_matches = points.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.intersects(poly)]

    # count number of firms within polygon
    x=pd.DataFrame({'grid_index': i,
                    'geometry':[Polygon(poly.envelope)],
                    'n_firms': len(precise_matches)})

    for poi in range(1,9):
        x['sic_'+str(poi)]=len(precise_matches[precise_matches.sic1==poi])
    return x

But I'm stuck in the next step which I think should be something like this

pool = mp.Pool(mp.cpu_count())

results = pool.map(grid_inter, [poly for poly in grid])

pool.close()

How can I make ithis work?


texts: I am trying to make a map using ArcMap 10.3 and adding Excel coordinates from a map. I originally converted addresses to lat/long using an online site. I have tried every suggestion I could find online and I can't figure how to get the points to show up. I've been able to upload the Excel file and add xy coordinates, but nothing shows up when I do. I only have 12 addresses/coordinate pair to add so I also tried to add them individually and couldn't figure out how. How can I fix this?


texts: Control to select data source on Leaflet map

texts:
I am building a map using leaflet.

My data source looks something like this:

var districts =
{"type": "FeatureCollection", "features":
[

{"type":"Feature",
"geometry":
{"type":"MultiPolygon","coordinates":[[[ [-75.228404,40.084708],[-75.2286,40.084875],[-75.230501,40.086024]//...
]]]
},
"properties":{"dist_num":"8", "gender":"M", "party":"democrat"}},

{"type":"Feature",
"geometry":
{"type":"MultiPolygon","coordinates":[[[ [-75.227985,40.084356],[-75.2281,40.0843],[-75.228241,40.084198],//...
]]]
},
"properties":{"dist_num":"5", "gender":"F", "party":"republican"}}, //....

];

I've got a working map that shows districts by party (democrat/republican) using this javascript:

var map = L.map('map').setView([39.955, -75.1647], 11);

L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: 'Map data &copy; OpenStreetMap contributors, CC-BY-SA, Imagery © Mapbox',
    maxZoom: 18
}).addTo(map);

var geojson;

function getColor(d) {
    return d == 'Democrat' ? 'blue' :
           d == 'Republican'  ? 'red':
              '#FFEDA0';
}

function style(feature) {
    return {
        fillColor: getColor(feature.properties.party),
        weight: 2,
        opacity: 1,
        color: 'white',
        dashArray: '3',
        fillOpacity: 0.3
    };
}

L.geoJson(districts, {style: style}).addTo(map);

Now I want to add a control where the user can switch from the fillColor in the function style from political party to display by gender. Does Leaflet have a built in method to do this? Or do I need to build my own in HTML/Javascript/Jquery?


texts: Different values when trying to use window in rasterio

texts:
I've been having a problem using windows in rasterio, which is probably a misunderstanding on my part.
My goal is to dump a specific data from a band into a csv file, first I tried to use a window but the results differ from the results that I see on QGIS (the results appear to be shifted 1 position).
So I decided to try to read the whole band, instead to use a window, and the result match with QGIS.
I think that I am doing something wrong when I create/use, what I tried to do was:

create a window from the bounds of the desired area:
bandWin = rasterio.windows.from_bounds(*boundingBox,band.transform)

get the transform matrix of that window:
win_transform = band.window_transform(bandWin)

Obtain the window data:
win_data = band.read(1, window = bandWin)

Get the desired window indexes:
indexes = rasterio.transform.rowcol(win_transform, [Xs], [Ys])

Get the values that I want:
csvData.append(win_data[indexes[0],indexes[1]][0])

The only differences between when I don't use a window is that instead of using win_transform I use band.transform and read the whole band.


texts: Prioritising values in accordance to various Columns

texts:
I previously asked for help in Prioritizing the Attribute table of a self made polygon layer
Prioritising polygons in QGIS

Now I want to change this prioritisation in accordance to multiple values of multiple columns.

I am looking for a nice way to add the logic

If "Area Size" > 10 ha then Prio = 1

If "Sun intensity" = high then Prio -1, otherwise Prio stays the same

If "Groundpoints" > X, then Prio +1, otherwise Prio -1

If "Restriction" = True, then Prio = 0

the minimum Priority possible is 1, because 0 determines the site of interest to be deleted and is only possible to be acchieved if "Restriction" = true.

The highest Priority possible is 3, which can either be determined through Area size (previous post) or additional factors (current post)

I am aware that I can achieve this command through endless WHEN, AND or OR Cases, however, I am looking for a clean and easy to understand solution.


texts: How do you classify unclassified pixels in a land cover geotiff tile?

texts: While printing out the raw pixel values of a GEOTIFF land cover tile I found 90 % of pixels have normal values i.e. those that the documentation describes (8 bit integers describing a land cover type). However several raw pixel values have been found to be zero. From this link - http://eoedu.belspo.be/en/guide/postclas.asp?section=3.8 I presume this means that they are unclassified.Has anyone performed these post-classification techniques such as modal filtering  using their own code or other software so that these unclassified pixels do get 'attached' to some land cover type ? If necessary I can upload an image of the necessary GEOTIFF file.


texts:
Is it possible to render a background color for text labels in geoserver? My current style is this:

  <TextSymbolizer>
    <Label>
      <ogc:PropertyName>fahrznr</ogc:PropertyName>
    </Label>
    <Font>
      <CssParameter name="font-family">Arial</CssParameter>
      <CssParameter name="font-style">Normal</CssParameter>
      <CssParameter name="font-size">12</CssParameter>
    </Font>
    <LabelPlacement>
      <PointPlacement>
        <AnchorPoint>
          <AnchorPointX>0.5</AnchorPointX>
          <AnchorPointY>-1</AnchorPointY>
        </AnchorPoint>
      </PointPlacement>
    </LabelPlacement>
  </TextSymbolizer>

Which produces the following text (with transparent background):

But what I want is the following (with green background):

Is this somehow possible with geoserver style? I tried

<CssParameter name="fill">#00FF00</CssParameter>
<CssParameter name="fill-opacity">0.7</CssParameter>

And

<CssParameter name="background">#00FF00</CssParameter>

Both do not work.


texts:
How to write a query in PostGIS  to check and fix case in picture:
In the element_polygons I store the geometry geom.


texts:
I am trying to write georeferencing information to a '.las' file using liblas and c++.

Here is my code:

   ofstream ofs;
   ofs.open("test.las", ios::out | ios::binary);
   liblas::Header header;
   liblas::Writer writer(ofs, header);

   liblas::Point point(&header);
   point.SetCoordinates(1.0, 2.0, 3.0);

   writer.WritePoint(point);
   header.SetPointRecordsCount(1);

   liblas::SpatialReference srs;
   srs.SetFromUserInput("EPSG:4326");
   header.SetSRS(srs);
   writer.SetHeader(header);
   writer.WriteHeader();

But, when I check the file using lasinfo tool I get:

 Spatial Reference:           None

If I check proj4, the information is stored, so I don't know why lasinfo doesn't print it:

Code  : std::cout << srs.GetProj4() << std::endl;
Output: +proj=longlat +datum=WGS84 +no_defs

However, when I use lasinfo tool with others las files not written by me I can see the WKT and geotiff output.

What am I doing wrong?


texts:
I would like to update the Map_Feature_ID field with a String " WP_" and then follow the "WP_" with the OBJECTID field.

So I will end up with "WP_1", "WP_2" all the way to the end in the Map Feature ID field.

I am not sure if there is a way to do it with a python script or should I use the field calculator instead?


texts: How to add diagrams to a layer in a PyQGIS standalone application?

texts:
In my standalone application, I show different layers based on queries into a canvas.
Now I want to add a layer with charts over by the values of an attribute.
Is it possible in a standalone application?
I tried "QgsDiagramRenderer" but without Results.
What should I do?


texts: QGIS - QML/HTML-Widget access Data from other layer

texts:
Is it possible to access data with a qml/html-widget of a other layer?
I want to show specific fields, which are stored in a second layer in some kind of table. I'm unsure where or how to begin.
If it's possible what widget-type would be better suited?
EDIT:
I have a geopackage with one point layer and a non-geom layer (just non-spatial data). The non-geom layer is a history of the point-layer. So when I update the data of the point layer, a trigger defined in the point-layer inserts automatically a history-entry with the former actual data of the point-layer.
Now I want to see in the point-layer (on a second tab) the history of the points in a table. In the history-table is the ID of the point-layer for the reference.
So I want to know, how can I access the history-layer from within an qml/html-widget of the point-layer.


texts: Calculating corners coordinate from World file(.jgw)


texts: Least cost path with one-off cost for highest elevation

texts:
I'm trying to write a script to calculate the least cost path between two cells on a raster. I'm reasonably familiar with the approaches used by tools such as GRASS GIS's r.cost and ArcGIS's CostDistance tools. However, there is one constraint I can't figure out how to include:

The path represents the route of a pipeline carrying water. An increase in elevation incurs a cost (due to the pumping required). However, for a given elevation this is a one-off cost. i.e. the pumping cost is only related to the height of the source and the single highest elevation the path travels through. If the path passes through a cell at 5m elevation, it doesn't cost anything to travel through any other cell at the elevation.

How could I go about incorporating this constraint in my analysis? I'm open to any solution, but preferably one that could be scripted (and extra points for free/open source!).


texts: How to reference to my in QT designer created buttons in python?

texts:
I have created an own QDialog (using QT Designer) with three push buttons. The names are pButton1, pButton2 and pButton3. After saving the file I created the ui file by using:

pyuic4 ui_xy.ui > ui_xy.py

I can see my created GUI in QGIS, but I have no idea how to get a reference to the buttons. When I use QMessagebox it is very easy to get to know wich button was pushed, but using the same way for this GUI doesn´t work.

Using QMessagebox the following code works:

ret = QMessageBox(..., QMessageBox.Ok)
if (ret = QMessageBox.Ok):
   DO SOMETHING
else:
   DO SOMETHING

When I use the standard GUI with the Ok and Cancel button, the following code works:

dlg = xyDialog()
dlg.show
result = dlg.exec_()
if result = 0:
...

How can I obtain the same result using my buttons? Replacing the 0 with pbutton1 (or 2 or 3) doesn´t seem to work.

Do you know a good solution or tutorial?

Thanks!

Best,


texts: Constructing polygons from GDAL Polygonize?

texts:
I have a geotiff file, which, for each cell ( pixel), it contains a discrete value ( you can think of this value as the "fitness" value of the cell). For the neighboring cells, they can consist of the same discrete value, but not always.

I want to create a polygon for all the adjacent cells that share the same discrete value. A polygon is defined by the following usual class ( in pseudo code)

class Point
{
    double X;
    double Y;
}

class Polygon
{
    List<Point> Pts;
}

The reason for this is that, for my rendering engine, I would rather render a few large polygons, rather than all the individual raster cells one by one for performance consideration.

I believe that GDALPolygonize is the way to go. However I don't quite know how to turn the output from the function into a list of Polygons that I can use.

The example from the cookbook straightaway makes use of a library to render the output file, which is not suitable for my purpose.

I mainly working with GDAL ( C++) or GDAL.Net (C#), but examples in other languages are acceptable as the syntax is quite close with C++ version.

Or is there other, more direct method that I can use, both in or out of GDAL?


texts: Convert over 360 degree range of coordinates to +/-180 longitude range in R?

texts:
I'm working with data on sea surface currents from ftp://podaac-ftp.jpl.nasa.gov/allData/oscar/preview/L4/oscar_third_deg The data has the latitudes and longitudes saved in the following format:

latitude  Size:481
units: degrees-north
long_name: Latitude
longitude  Size:1201
units: degrees-east
long_name: Longitude

head(lon)
[1] 20.00000 20.33333 20.66667 21.00000 21.33333 21.66667
tail(lon)
[1] 418.3333 418.6667 419.0000 419.3333 419.6667 420.0000

How can I convert this into latitude-longitude units in R?

The lon object goes from 20 to 420 whereas I want it to be from -180 to 180.

I have seen Convert North and East coordinates to longitude and latitude in R at Stack Overflow but I am unable to adapt it to my setting.

x <- nc_open("oscar_vel1992.nc")
# lon
lon <- ncvar_get(x,"longitude")
#lat
lat <- ncvar_get(x,"latitude")
z <-as.numeric(char2dms(lat))


texts: How to copy all files with extension type to a folder?

texts:
What is the best way to go about this? I'm able to do it with a single file as below, but what if I would like to copy all files with .txt to a folder, and all files with .docx to another folder. I tried changing in_data to just .txt but this did not work.

# Import system modules
import arcpy

# Set workspace
arcpy.env.workspace = "C:\Users\X220\Documents"

# Set local variables
in_data =  "test.txt"
out_data = "C:\Users\X220\Documents\output\Copy.docx"

# Execute Copy
arcpy.Copy_management(in_data, out_data)


texts: How to generate route to pass through polygons?

texts:
A twist on a standard problem. I'm looking to plot a route through all the wards in London (649). I have all the wards as a shapefile and ideally i wanted to let the computer generate a spiral route staring at the outer edge of london working inward to the centre. My thoughts of how to approach;

Assign each polygon a centroid.
Generate route to pass though points, outputting a line
Create a buffer round that line (to select all nearby minor roads)
Generate final route by rerunning algorithm along the roads

Does that approach seem valid with the tool sets available?

Image of polygons as below


texts:
this is a novice question

I have an attribute field that's missing a piece of information, i'd like to amend the attribute field like below:

This is the current attribute field:

"V:\Building & Development...."

and i'd like to add "\Engineering" to attribute field

"V:\Engineering\Building & Development...."

I have to make this change to 1441 attribute fields, not sure how to write the python script in the field calculator.  I'm doing this on ArcGIS 10.2 for Desktop


texts: default arcgis account credentials for 10.2 AMI

texts:
I'm installing the GeoEvent Processor on an ArcGIS 10.2 AMI Server and the install dialog is asking for the password for the arcgis account.  Because the server was created with Cloud Builder, I never set a password for the account.

Is there a default password for the arcgis account?


texts:
I can't seem to get the mosaic created for this.

 import arcpy

 arcpy.env.workspace = r"C:\projectTemp\FloodAnalysis\output\2yr\raster"
 arcpy.env.overwriteOutput = True
 srName = r"C:\projectTemp\FloodAnalysis\output\2yr\raster"
 outPut = r"C:\projectTemp\FloodAnalysis\output\2yr"
 rasterName = "2ft.img"
 bands = 1
 method = "First"
 clMap = "First"

 arcpy.MosaicToNewRaster_management(srName, outPut, rasterName, "", "", "", bands,      "First", "First")


texts: Finding longest straight line within polygon in QGIS

texts:
I'm looking for an algorithm to draw max distance straight lines contained in a polygon (without crossing polygon boundaries). How can I perform it?
As you can see some polygons are very oddly shaped.

The area I am working with is millions of acres and contains hundreds of thousands of individual polygons. From what I can tell the polygons are pretty clean (no holes, overlaps, etc.)
The practical use of this is to filter what bodies of water are capable of having a floatplane land and take off on them.
The only criteria I have for a solution is that I can perform it in QGIS and that the lines are straight.


texts: Reading values from raster, and creating column of those values in Attribute Table of shapefile?

texts:
I have a map of polygons from a shapefile which lost its numbered labels in the conversion process from .pdf to .shp. In order to make these labels visible again, I creased a raster from a .png file, and georeferenced this layer onto its identical polygon layer.

My goal is to create a column of values in the Attribute Table of the polygon layer that includes the label from the raster layer for each specific polygon. Is there a simple way to auto-fill this column? Although the labels are numbers, the labels can be read as text to prevent any formatting error that could arise.

I need the polygons to have the labels from the raster layer because I have information (land-use, price, etc.) for each polygon based on these labels.

In the attached screenshot below, the selected polygon has ID "2917" and should have its CTS column filled with the value/string "1572/8".


texts: How can I Create Line Feature using two separate point features in MapInfo

texts:
I have two MapInfo Tab files as Point feature. both the layer have the same unique ID. what I was trying to do is to connect both point using the same ID  as a Line feature. so the end result will be "Line Feature" created by two separate point  tab file. please help me to resolve this issue using MapBasic. currently I am using MapInfo 15.2.

Above image states two separate point feature that I was tying to connect which have same id field.


texts: ST_Intersection preserving overlapping edges


texts: Uber H3 DGGS: Table of largest hexagons per resolution?

texts:
I'm evaluating the use of H3 as a DGGS for a project, and so far I've been unable to find any reference of the area of the largest and smallest hexagons in the grid at each resolution. The "official" resolution table provided on the H3 website lists only the average edge lengths and areas, and I know that the maximum area of a hexagon in a given resolution is approximately 1.6 times the area of the smallest, but I need precise values.

Is anyone aware of a table that lists these, or is this likely something I'll have to do myself by just generating grids at each resolution and finding the maximum and minimum areas?


texts: Taking data from a CSV file to create points on a shapefile

texts:
I am trying to import data from a CSV file and write a script that will generate points on a new shapefile.  I am running into issues for when I try to update the InsertCursor. I get an error of "not enough quota is available to process this command".  I am fairly new to Python; does anyone have any idea of what is wrong with my code?

#Import a csv file and take the coordinates of each ASOS and create a polypoint output
import csv
import arcpy
import os
from arcpy import env

#Lets you write over the shape file if needed
arcpy.env.overwriteOutput = True

targetFolder = "C:\\GEOG485\\FinalProject\\OutPut"
env.workspace = targetFolder

#csv file that we are importing
asosFile = "C:\\GEOG485\\FinalProject\\isd-history.csv"

targetState = 'NJ'

try:
    #Create the featureclass ASOS shapefile, 4326 is the factory code for the required coordinate system of the file
    asosStations = arcpy.CreateFeatureclass_management("C:\\GEOG485\\FinalProject\\","ASOS","POINT",'','',"ENABLED",4326)
    print ("success")
    #add a new field called "Station Name", "CTRY", "STATE" and "END"
    arcpy.AddField_management(asosStations,"STATION_NAME", "TEXT")
    arcpy.AddField_management(asosStations,"CTRY","TEXT")
    arcpy.AddField_management(asosStations,"STATE","TEXT")
    arcpy.AddField_management(asosStations,"ICAO","TEXT")
    arcpy.AddField_management(asosStations,"END","TEXT")
    print ("success")

except:
    print ("error")

#asos dictinary
asos = {}

try:
    #Use the with statement to open the csv file and r to read mode
    with open(asosFile, 'r') as csvv:
        #Reading the csv file with DictReader which takes the first row of the file and uses them as keys
        reader = csv.DictReader(csvv)
        #Loops through the csv file
        for row in reader:
            #Dictionary to see what asos's have a key in the dictionary and append the asos location with the coordinates
            if row['ICAO'] not in asos:
                asos[row['ICAO']] = []
            asos[row['ICAO']].append((row['LAT'],row['LON']))
            print ("success2")
except:
    print ("error2")

try:
    #Takes the points that we just appended to each asos station and adds them to the shapefile
    with arcpy.da.InsertCursor(asosStations, ["ICAO","SHAPE@XY"]) as cur:
        #Loop checks for keys (file headers) and values (asos's inside the dictionary)
        for key,value in asos.items():
            print key,value
            #creates a point with the points for each asos station in the shape file
            for coords in value:
                point = arcpy.Point(*coords)
                #Takes all the keys with the same names and adds all the associated points to them
                row2 = (key,point)
                cur.insertRow(row2)
                print ("success3")
except:
    print ("error3")


texts: Change Leaflet-Control-Window modal content by clicking geojson feature

texts:
I have a leaflet map using Leaflet-Control-Window to produce a modal pop-up with dynamic content.

Using their documentation I was able to produce a popup when clicking on features of my geojson file, but I cannot figure out how to change the content of that modal based on my geojson file.

This calls the modal:

var geojson = L.geoJson(json_Locations3, {
  style: style,
  onEachFeature: onEachFeature

}).addTo(map)

    geojson.on('click',function(geojson){

      var win =  L.control.window(map,{title:'Hello world!',maxWidth:400,modal: true})
              .content('Info:' + Feature.properties.Comments)
              .prompt({callback:function(){alert('This is called after OK click!')}})
              .show()
});

This is a part of my geojson file:

var json_Locations3={
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
{ "type": "Feature", "properties": { "id": 2001.000000, "Comments": "Laguardia Houses. Filler text and description.", "Address": null }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ -73.98707531157136, 40.71279019297274 ], [ -73.990326197834079, 40.710554740197978 ] ] ] } },
{ "type": "Feature", "properties": { "id": 1002.000000, "Comments": "Essex Crossing", "Address": "236 BROOME STREET" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -73.988559744585103, 40.718138957271904 ], [ -73.988629841437202, 40.718003100250257 ], [ -73.988805496163607, 40.718056044119535 ], [ -73.988939514855289, 40.71779629683828 ], [ -73.98937138176332, 40.717926433649048 ], [ -73.98917121142054, 40.718314400521201 ], [ -73.988867202969402, 40.718222771755336 ], [ -73.988863247936536, 40.718230436389028 ], [ -73.988559744585103, 40.718138957271904 ] ] ] } },
{ "type": "Feature", "properties": { "id": 1003.000000, "Comments": "Essex Crossing Test Description", "Address": "145 CLINTON STREET" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -73.986825413362538, 40.717732897085796 ], [ -73.987087772457556, 40.717230210278437 ], [ -73.987759556669246, 40.717435703706677 ], [ -73.987494416232039, 40.717949651831589 ], [ -73.986819055387571, 40.717745220879294 ], [ -73.986825413362538, 40.717732897085796 ] ] ] } },

Any suggestions? I am really a beginner with javascript and struggling with different solutions, but not having luck.


texts: Geoserver Clustering. where can i find activemqBroker.war?

texts:
I am completely new to clustering setup and I am following this tutorial to setup geoserver cluster. But I am stuck at

  "HOWTO configure ActiveMQ broker"

step, problem is i cannot find this activemqBroker.war Its mentioned here

  "Deploy the produced activemqBroker.war in your tomcat instance "

Where can i find this activemqBroker.war ?


texts: Why does 3d Bounding Box &&& check return true when I use ST_3DMakeBox?

texts:
I'm just wondering why this &&& is returning true when I use ST_3DMakeBox?

SELECT ST_MakePoint(0,0,0) &&& ST_3DMakeBox(
  ST_MakePoint(-1,-1,-1),
  ST_MakePoint(1,1,-1)
);

When these &&& returns false

SELECT ST_MakePoint(0,0,0) &&& ST_3DMakeBox(
  ST_MakePoint(-1,-1,-1),
  ST_MakePoint(1,-1,-1)
);

SELECT ST_MakePoint(0,0,0) &&& ST_3DMakeBox(
  ST_MakePoint(-1,-1,-1),
  ST_MakePoint(-1,1,-1)
);

Going through the X-Y, I can see a non-intersection in those ranges makes the 3d-test return false. However, a non-intersection in Z, even with &&& is still returning true.

    postgis_full_version
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 POSTGIS="2.2.2 r14797" GEOS="3.5.1-CAPI-1.9.1 r4246" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 2.1.2, released 2016/10/24" LIBXML="2.9.3" LIBJSON="0.12.1" RASTER
(1 row)

Just an Ubuntu 17.04 stock install.


texts: Strange behaviour viewing world map with EPSG:4326

texts:
I'm trying to visualise/convert EPSG:3857 map tiles to EPSG:4326, but am experiencing strange behavior in that I don't think the most zoomed out view of the world is using the correct CRS/Auth ID.

EDIT: this behavior only seems to occur when using XYZ map tiles.  If testing with the Natural Earth vector data, it is displayed correctly at all zoom levels.  Overlaying the XYZ map tiles over the Natural Earth data show the same disparity (i.e. the tiles seem top-bottom squashed).

To replicate this (strange to me) behavior:

Create a new QGIS 3 project
Load in XYZ tiles (e.g. one of: ESRI Standard, Google Maps, OpenStreetMap standard)
Set Project CRS (in bottom right) to EPSG:4326
(Optional) Right click raster tile layer, and 'Zoom to Layer'
If you then zoom in (with a single scroll wheel click) to say Brazil, or Africa, its shape significantly changes.  Subsequent zooms seem to work as expected.
Also, plotting EPSG:4326 data at the most zoomed out level seems to put things in the wrong place. Upon zooming in, they jump to the right place.

My questions:

Why is the CRS/Auth ID of the most zoomed out image (Fig 1) not EPSG:4326?
Is this behavior to be expected? Is it a bug in QGIS? Is there a workaround (to make the initial overall world view correct)?

Fig 1: Example of zoomed out image:

Fig 2: Example after single scroll wheel click, to zoom in:


texts: Creating QGIS legend showing different size points from set scale field?

texts:
I have created a map with different sized points based on another attribute using the set scale field in QGIS v2.8.4. See map image.
I am trying to create a legend like this previous thread showing different size bubbles with the range of attributes. How can I do this?
Resizing symbols in QGIS Print Composer legend?


texts: I regularly use land registry data from https://www.gov.uk/government/collections/download-inspire-index-polygons however at present the .gml, when adding a vector layer, renders points and not polygons (as it usually does). Could anyone advise how this could be rectified?


texts:
I need some basic instructions on how to create a network dataset in ArcMap 10.1. I do realize there are esri tutorials out there but I'm having difficulty understanding everything. I'm just looking for a little help to get started.

What I want to do: perform Network Analysis on a potential new hospital site - I want to see how this new site will compare to an existing hospital location in terms of how many people it covers, accessibility, etc.

So first I need to create a Network dataset. I currently have streets and public transit route shapefiles, and I created 2 polygon features for the old and the new hospital location. What else do I need? Where do I go from here?


texts: Find mean / median point based on ID attribute and point density (hot spots)

texts:
I have a point feature class within ArcGIS Pro, with 10.000 features, data are structured like this:
ObjectID |  MemberID | Lat | Long
1             ID_1     Y1    X1
2             ID_1     Y2    X2
3             ID_1     Y3    X3
4             ID_2     Y4    X4
5             ID_2     Y5    X5
...
10.000        ID_n     Yn    Xn

Based on:

MemberID
Density of points (sometime I have points very far each other and I want to give a low weight to the outliers)

I would like to have as output:

A point feature class, in which each point represents the mean/median/avg/bestRepresentation of location for each MemberID.

I already tried with the centroid, but I don't like the result as I would also like to take into consideration the density of points (and give a low weight to the outliers)
Is there a specific tool to do this or should I built a custom one?
Once I receive suggestions I could start to write some code below.


texts: Licensed File Geodatabase to unlicensed Geodatabase?

texts:
Can anyone help me how to convert a licensed File Geodatabase "Export Deny" encrypted in ArcGIS to "allow export" if i have only File Geodatabase and  file without license definition  file.

I am in great trouble because maximum feature is disabled in the licensed File Geodatabase.


texts: What is the difference between Key-Point, Tie-Point in Photogrammetry?

texts:
I'm looking into the reports & the Mechanism of UAV Photogrammetry software like Pix4D and AgiSoft.

These two words are mentioned in a lot of support documentation, but they are defined no-where.

From my studies, I remember that tie-points are points marking the same physical feature in two or more images.

What does key-point mean? Is it the same?


texts:
Just ran my first iterations of data-driven pages. It was more successful than I planned, but I'm running into some issues.

Main problem is rotating the page to the extent of the polygon. There are two areas where this is a large problem because it minimizes the size of otherwise important segments.

Here is my DDP setup:

I calculated the orientation of the polygons with the Grid Conversion to Angle tool:

Here is a section of the index sheets along the highway with the resulting orientation (and a x100 attribute of orientation I created called Rotation, which still doesn't appear to be a large enough value for a degree angle):

However, this is how the resulting DDP exported page (with index page turned on for comparison) turned out:

Everything is rotated to true north as opposed to the orientation of the index page. Additionally, it creates the output to the extent of the polygon vertices as opposed to rotating to match the polygon's general direction.

What am I doing wrong?


texts: Looking for WFS-T which data I can edit in QGIS for testing purposes

texts:
I'm in a situation where I have to make sure, that a WFS-T solution for qgis works on my office-desktop here in a secured network, before the responsible IT-department can creates an environment with webserver, database and qgis-server for that.

Has anyone an idea, where I can find a WFS-T which I can edit and save back to prevent reactions from IT like "make sure it works in this network with that proxy before wasting my time"?


texts: Creating Network Dataset programatically

texts:
I'm trying to code to create a network dataset using a single line feature class with minimal attributes
I'm looking for a nice tutorial which can guide me through the process.I have already looked at the examples provided by esri (http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//000100000466000000) but they are really complex and specific to the question, its really hard to get the concept from it.


texts: Importing GeoJSON without duplicating values in properties

texts:
Here is the issue I am having, it should be fully reproducible. I am attempting to import a GeoJSON into ArcGIS Pro. The issue is that if one of the features has a property that the other doesn't, instead of it being null in the attribute table, ArcGIS uses the last features value as the new features value.
I included a tiny GeoJSON that should reproduce the issue. I am importing it through json -> features.
{
    "features": [{
        "geometry": {
            "coordinates": [-74.2462, 43.4674],
            "type": "Point"
        },
        "properties": {
            "area": "",
            "height": "25",
            "id": "13784",
            "link": "https://digthefalls.com/auger-falls-wells/",
            "name": "Auger Falls",
            "node_type": "waterfall",
            "waterway": "Sacandaga River"
        },
        "type": "Feature"
    }, {
        "geometry": {
            "coordinates": [-74.6567, 43.6012],
            "type": "Point"
        },
        "properties": {
            "area": "West Canada Lake Wilderness",
            "id": "13785",
            "name": "Brooktrout Lake Shelter",
            "node_type": "lean-to"
        },
        "type": "Feature"
    }],
    "type": "FeatureCollection"
}

When imported, the second feature in the GeoJSON will have the same link as the first feature in the GeoJSON. I would like the second features link to be null.


texts: I am looking for a list of mapping libraries?

texts:
I am looking for the top APIs that I can use to display maps from a number of map servers/sources. (like openlayers and openscales). So are there any other popular APIs for other languages than javascript and flex (C++,Java,C#... ) ?
(it doesn't matter if it is free or not)


texts: filter out atlas pages with no data in QGIS 3

texts:
I am printing maps using the atlas to display layers with point data.  Some of the areas have no points in them.

Is it possible to suppress these using the filter options.  I could not figure out how.

i.e. don't generate the map unless the area contains one or more points in the area.


texts: Can't open world.qgs as in Training Manual 6.1.1

texts:
I’m working through the training manual located at https://docs.qgis.org/3.16/en/docs/training_manual/index.html. I have QGIS 3.16 running on OS 10.15.7 (Catalina). The following describes my problem:

Module: Vector Analysis 6. 1.1 Follow Along: Projections.

The instructions say to save the current map and then open a map of the world under exercise_data/world/world.qgs.  I opened a new project. Then I selected Open Data Source Manager and, Since Module 6 pertains to Vector Analysis, I selected Vector on the left toolbar, and then navigated to the world.qgs file. A red banner appeared with the message that this was not a valid data file. I tried opening it without opening a new project with the same result. I then opened the previous project and tried open the file in that project. I also tried opening the file as raster data and as mesh data. The result was the same. I would like to know where I have gone wrong.


texts: Is it possible to carry out vector editing functions in QGIS? In this case I have an imported shapefile of contour lines. I wish to be able to create a copy of a line, in-situ, trim it, then create a number of offsets at a specific spacing. In short, basic CAD tasks.


texts:
I have a .shp file, where I want to connect all points (approx. 100) by lines.
I have a data field with the order (i.e. 1, 2, 3, ...).

I have tried PointConnector and Points2One, maybe I am using them wrong, but I am not able to get any results.

In my searches, I've found mention of the plugin points to paths, but that seems to be gone?

I'm using QGIS 2.6.1


texts:
I would like to calculate coefficients of a 2D affine conversion between two images.
Let me know of online pages that can help make programs using C or C++ or provide an application. I need a program or application which can calculate coefficients of an affine conversion using coordinates of corresponding points.


texts: Proper aliasing of spatial columns in SQL

texts:
Referencing the alias length_road in the WHERE clause does not seem to be allowed. Does this query as written cause the ST_LENGTH operation to be run more than once per record, and if so, how can that be avoided?

SELECT
  roads_nat.*,
  ST_LENGTH(roads_nat.road_geom) AS length_road
FROM
  `bigquery-public-data.geo_us_roads.us_national_roads` AS roads_nat
WHERE
  ST_LENGTH(roads_nat.road_geom) > 1000 -- vs. length_road > 1000


texts: tmap max categories being weirdly combined

texts:
I am trying to plot homeless numbers in the UK using tmap. However, I am coming into some problems with the max.categories for tmap_options.

Data

EnglandWales <- geojson_read("https://opendata.arcgis.com/datasets/8edafbe3276d4b56aec60991cbddda50_2.geojson", what = "sp")

Homeless data:
https://data.cambridgeshireinsight.org.uk/sites/default/files/P1E-national-homelessness-CLG-tab784-to-201718-csv.csv

So I merge the data:

EnglandChildDeprivationDataMap<-merge(EnglandMap,
                                      EnglandHomelessData,
                      by.x="lad15cd",
                      by.y="ONS code",
                      no.dups = TRUE)

and then I plot using:

tm_shape(EnglandChildDeprivationDataMap) +
  tm_polygons("2017/18: Numbers accepted as homeless and in priority need total",
              main.title = "Homelessness Count In The Uk",
              style="cont",
              palette="Purples",
              midpoint=NA,
              title="2017/2018 Recorded Number Of Homeless People") +
  tmap_options(max.categories = 15)

From this, I get the warning message:

"Number of levels of the variable "2017/18: Numbers accepted as homeless and in priority need total" is 199, which is larger than max.categories (which is 15), so levels are combined. Set tmap_options(max.categories = 199) in the layer function to show all levels"

and my categories get merged like such:

Why are my categories being merged like this, and how can I make it more understandable?


texts: Performing water valve isolation trace using ArcPy?

texts:
The utility that I am working with for a couple of weeks has access to the ArcFM Water Valve Isolation tool which enables them to point from ArcMap at a location on their SDE Geometric Network and be shown the water valves with NORMALPOSITION OPEN that would need to be CLOSED to isolate that part of the network in the event of a pipe burst.

I am not strong on ArcFM, but I what I am hoping to do from ArcPy (with ArcFM licensing available) is to take a list of 2500 locations on the geometric network, and have an isolation trace return a list of all valves that any of these identify as being used as potential isolation valves from bursts at those locations.

Doing 2500 isolation traces from ArcPy should be easy if I can do one isolation trace from ArcPy.

Is it possible to perform an isolation trace from ArcPy with ArcGIS 10.2.1.01 (Standard level license; with or without ArcFM)?

From @Midavalo's answer to Finding valves to turn off for water line using ArcMap? I was hoping that ArcGIS for Water Utilities download might offer some suitable Python code but when I downloaded it the Valve Isolation tool seems to use ArcObjects instead.

I have also been looking at the Trace Geometric Network tool which is ArcPy accessible but does not offer the option of an isolation trace.


texts: Checking whether Raster has Attribute Table using ArcPy?

texts:
I can't seem to find a Describe property or any tool that will tell me whether a raster has an attribute table.

I can do it with a simple try/except statement on a test cursor or a table-specific tool, but those return relatively generic RuntimeError and ExecuteError, respectively. That seems messy and insufficient.

Is there a property/tool that will give me a simple boolean indicating whether the raster has an attribute table?

It's not that I want to ensure that there is one; it's that I simply want to test whether there is one.


texts: Which feature has highest value within extent

texts:
I have a layer in which I need to find out the maximum value of a given field within the current map extent, and then return back the value from a different field.

e.g. Find state with highest population within current extent, and return the state name.

This is a hosted service on ArcGIS Online, so I thought the right way of going about this was to make use of outStatistics to determine the maximum value of my field using the max property.  However, even if I specify the outFields to include the field I need to report back on, all I get back is the result of the outStatistics.

For example, see this GET request.

Out Fields: POP2007,STATE_NAME,NAME
Extent: <an extent object>
Output Statistics: [{"statisticType":"max","onStatisticField":"POP2007","outStatisticFieldName":"maxPop"}]

Result:  maxPop: 1483460

Desired response:

Result:  maxPop: 1483460, POP2007: 1483460, STATE_NAME: Middlesex

What is the most efficient/elegant way of calculating the maximum value of one field within the current extent, finding that record, and returning a value from a different field?

In the above example, we can forget about extent, I just want it to return the State Name that has the maximum population value.  Please also disregard the fact that counties are involved.


texts: QGIS 2.14.4 certificate error after installation

texts:
When QGIS 2.14.4 is being deployed to a new PC, we get the following error message:

Unable To Get Local Issuer Certificate: The issuer certificate of a locally looked up certificate could not be found
Certificate Untrusted: The root CA certificate is not trusted for this purpose

How do we resolve this issue?


texts: Postgis geometry creation of polygon from geojson

texts:
I am working on postgis and was ran below query but i have my doubts

Is the below representation valid ?

As far as I know the polygon must have 4 coordinates.

SELECT  ST_GeometryType  ( ST_GeomFromGeoJSON('{
        "type": "Polygon",
        "coordinates": [
            [
                [
                    10.798460841178894,
                    27.7146573614866
                ],
                [
                    18.005492091178894,
                    23.112072196617056
                ]
            ]
        ]
    }'));

What can be work around for the validation ?

Can we throw error before insert in such case ?


texts:
I have the following data:

The pink represents digitized hiking trails done on a different software that is completely inaccurate compared to raster images. They do however have attributes with the hiking trail names.
The rest of the features represent true and accurate street and walking path features. I need to somehow make it so that the features similar to the pink ones are extracted and have the attributes of the digitzed hiking trails.

I am using ArcGIS 10.1 and have no Spatial Analyst and Network Analysis.

FME solutions are also welcome


texts: Wrong number format when creating SLD from YAML via REST API

texts:
In our workflow we are creating style definitions in YAML locally and then upload these styles via GeoServer REST API to the server. A sample YAML file looks like this:

name:  Test Style
title: test style description
abstract: Used for testing
feature-styles:
- rules:
 - title: raster
symbolizers:
  - raster:
      opacity: 1.0
      color-map:
        type: values
        entries:
        - ['#e20374', 1.0, 1, lorem ipsum]

The creation on the server works fine. However, when we try to apply this style to a layer, we get the following error:

2017-03-01 17:40:36,234 INFO [geoserver.wms] -
Request: getServiceInfo
2017-03-01 17:40:36,240 ERROR [geoserver.ows] -
java.lang.NumberFormatException: For input string:"java.awt.Color[r=226,g=3,b=116]"
at java.lang.NumberFormatException.forInputString(NumberFormatException.java:65)
at java.lang.Integer.parseInt(Integer.java:580)
at java.lang.Integer.valueOf(Integer.java:740)
at java.lang.Integer.decode(Integer.java:1197)
at java.awt.Color.decode(Color.java:729)
at org.geoserver.wms.legendgraphic.LegendUtils.color(LegendUtils.java:494)
at org.geoserver.wms.legendgraphic.LegendUtils.color(LegendUtils.java:514)

In GeoServer when I check the style directly in the GUI, this is what I see:

name:  Test Style
title: test style description
abstract: Used for testing
feature-styles:
 - rules:
  - title: raster
    scale: [min, max]
    symbolizers:
     - raster:
         opacity: 1.0
         color-map:
           type: values
           entries:
           - ['java.awt.Color[r=226,g=3,b=116]', 1.0, 1, lorem ipsum]

As can be seen in this example, the #e20374 color code is converted java.awt.Color[r=226,g=3,b=116], which of course can't be converted to a number.

Question:

What causes this wrong conversion?
What needs to be done to ensure that the correct number string is applied?


texts: Apply transparency on a WFS layer (PyQGIS)

texts:
I'm trying to make a wfs layer transparent, but it gives me a mistake, but I can't solve it until now.
This is my code:
uri ='http://gstore.unm.edu/apps/rgis/datasets/6ca5428a-a78c-4c82-8120-da70dc92f2cc/services/ogc/wfs?SERVICE=wfs&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=tl_2010_35_state10'

wfslayer = QgsVectorLayer(uri, "MEXICO", "WFS")

QgsProject.instance().addMapLayer(wfslayer)

So far everything is correct
 renderer = wfslayer.rendererV2()

In this function it gives me the mistake
I don't know if it's a problem with the version (QGIS 3.16.2) of the program or with the code.
It no longer lets me put the rest of the code.
symbol = renderer.symbol()
symbol.setAlpha(0.5)


texts: How to Raster Layer checked status change?

texts:
I make plugin with use QListWidget.
Items of QListWidget is Raster Layer name.
When selected items, same name layer checked status changing.

current my code is

items = self.dlg.listWidget.selectedItems()
 layers = QgsMapLayerRegistry.instance().mapLayers()
 for name,layer in layers.iteritems():
  for item in list(items):
   if unicode(layer.name()) == item.text():
    utils.iface.setActiveLayer(layer)
    utils.iface.mapCanvas().setCurrentLayer(layer)

but it is not apply.
What should I do?


texts: geom_vline and scale_x_discrete dont fit with each other

texts:
I am a newbie here.

1, The aim is to plot a graph about the mean NDVI value during a time period (8 dates were chosen from 2019-05 to 2019-10) of my study site (named RB1). And plot vertical lines to show the date with a grass cutting event.

2, Now I had calculated the NDVI value for these 8 chosen dates and made a CSV file.
(PS. the "cutting" means when the grassland on the study site has been cut, so the corresponding dates should be show as a vertical line, using geom_vline)

infor <- read_csv("plotting information.csv")
infor
# A tibble: 142 x 3
   date         NDVI cutting
   <date>      <dbl> <lgl>
 1 2019-05-12 NA     NA
 2 2019-05-13 NA     NA
 3 2019-05-14 NA     NA
 4 2019-05-15 NA     NA
 5 2019-05-16 NA     NA
 6 2019-05-17  0.787 TRUE
# ... with 132 more rows

3, the problem is, when I do the ggplot, first I want to keep the x-axis as the whole time period (2019-05 to 2019-10) but of course not show all dates in between, otherwise there will be way too much dates show on the x-axis). So, I do the scale_x_discrte(breaks=, labels=) to show the specific dates with NDVI values.

Second I also want to show the dates that the grasses were cut geom_vline.

BUT, it seems like the precondition for scale_x_discrte is to factor my date, while the precondition for geom_vline is to keep the date as nummeric.
 these two calls seems to be contradictory.

y1 <- ggplot(infor, aes(factor(date), NDVI, group = 1)) +
  geom_point() +
  geom_line(data=infor[!is.na(infor$NDVI),]) +
  scale_x_discrete(breaks = c("2019-05-17", "2019-06-18", "2019-06-26", "2019-06-28","2019-07-23","2019-07-28", "2019-08-27","2019-08-30", "2019-09-21"),
                   labels = c("0517","0618","0626","0628","0723","0728", "0827","0830","0921")))

y2 <- ggplot(infor, aes(date, NDVI, group = 1)) +
  geom_point() +
  geom_line(data=infor[!is.na(infor$NDVI),]))

when I add the geom_vline in the y1, vertical lines do not show on my plot:
y1 + geom_vline

when I add it in the y2, vertical lines were showed, but the dates (x axis) are weird (not show as the y1 because we donot run the scale_x_ here)
 y2 + geom_vline

   y1 +
      geom_vline(data=filter(infor,cutting == "TRUE"), aes(xintercept = as.numeric(date)), color = "red", linetype ="dashed")


texts: How to overlay polygon shapefiles with different attributes?

texts: I have 2 shapefiles of the entire state of PA. One shapefile has median income for each county. The other has number of flu cases for each county. I would like to know, how should I combine them? My ultimate goal is to clip select counties from a combined shapefile.