texts: Create a raster that measures distance of two polygons

texts:
I'm currently working on a project where I have 2 datasets/layers that comprises boundary polygons (vector) of farms. One has been digitized from community mapping (overlaying a UAV image) the other one is collected using ArcGIS Collector. Now what I want to do is to know the distance between these two layers. For example Farm A boundary using UAV vs Farm A boundary collected by ArcGIS.

My goal here is to provide summary statistics on how accurate and fast is it to delineate polygons using a geo-referenced image without having to go to the site itself and collect data.

I'm using ArcMap 10.5 and wondering what tool (ArcGIS toolbar) I can use for this? My co-worker suggested that I can create a raster that would measure the distance between the two. Then clip the raster and calculate summary statistics.

I'm not really sure what to do based from her suggestion


texts: Why does import arcview/arceditor give gp.setProduct() RuntimeError: ERROR 999999: Error executing?

texts:
I am developing a script that will be used in ArcGIS Desktop 10.2 of Standard (formerly called ArcEditor) license level and thus importing arceditor to make sure I don't use GP tools that require Advanced (formerly called ArcInfo) license level (which I have) while developing the script.

I get the ERROR 999999: Error executing function when running the script in PyScripter. As outlined in the help, one has to import arceditor before importing arcpy. When running this single line of code:

import arceditor

When executing this single line of code, I get the C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arceditor.py file open in PyScripter and it seems that it has troubles setting product (this line is highlighted):

gp.setProduct("ArcEditor")

I have tried re-opening PyScripter, same error. import arcview gives the same error.
Did I get it right that one can set using Standard (formerly called ArcEditor) license while having an ArcGIS Desktop licensed with Advanced (formerly called ArcInfo) installed on the machine? If not in this way, how then?


texts: Working with NaN values in channel network in QGIS?

texts:
I have downloaded the conditioned 15s HYDROSHEDS dem of asia and after clipping the desired extent proceeded to evaluate the catchment area. However the next step of channel extraction renders NaN values.

Why does NaN value occur in the conditioned DEM?

Can initiation threshold can lead to such errors?


texts: Display two side-by-side maps in Google Earth Engine

texts:
Is there any way I can display two maps side by side on Google Earth Engine?

The same way leaflets can have a side-by-side plugin. Trying to find the widget to do so but can't find info.


texts: Applying different impendance by Break field for facilities in service area network analysis?

texts:
When I try to use my Break_Miles field to set own multi-impedance for every single facilities, it seems no effect but only use the default breaks. How can I apply my Break_Miles values to network solving?

According to http://desktop.arcgis.com/en/arcmap/latest/extensions/network-analyst/service-area.htm#GUID-21ADEC62-F784-4180-8D4D-547FB50621FD, I set this field as text field, the value is looks like "1 3 5" or "3 5 10", numbers divided by space means breaks value for each facilities. The Layer Properties also recognized my field as Breaks.

Any mistake I made?


texts: Programs/ languages to process netcdf files "fast"

texts:
I am using R to loop through thousands of NetCDF files (downscaled climate data), and while I am able to do it using R, the code is too slow. I am looking for a way to speed it up and to maintain functionality.

What I need is a program or platform where I can (1) import netcdf data, (2) import a polygon layer, (3) get the mean for the polygon layer, and (4) export output to csv file. The following is a general code workflow:

shape <- readOGR(dsn,"shape") ## read in shapefile
files <- dir("in directory", recursive=TRUE, full.names=TRUE, pattern="\\.nc$") ### get file paths to netCDF files
out <-NULL
for (i in 1:length(files)){
brick.tmp <- files[i]
    for (j in 1:nlayers(brick.tmp)){
    val = extract(brick.tmp[[j]],shape, fun = mean)
    out <- rbind(out,val)
}
}

I would prefer to stay in the R environment, but I'm not sure how to speed up the code.

I have attempted to use several programs and have looked through the list here:
https://www.unidata.ucar.edu/software/netcdf/docs/software.html

However, I thought I would leverage everyone's collective knowledge and see if I can narrow down the potential choices.


texts: Removing slivers from IPolygon using ArcObjects?

texts:
ArcGIS Desktop Basic 10.1 with Spatial Analyst and 3D Analyst

I have multiple neighboring polygons which I unionize. On occasion the result has slivers where the multiple polygons meet. I have been searching on the simplest way to get rid of these slivers but it seems like all tools work on features or feature classes. The result is a IPolygon element which I am drawing to the screen so don’t use the IFeature or IFeatureClass objects.

The before and after looks like this:

The end result should be one complete polygon with no slivers. The dimensions of the polygons are in 2 to 3 miles and the slivers are very tiny. I've zoomed in very close so estimate the slivers to be about ~5FT which is negligible for my purposes.

To unionize I have tried:

ITopologicalOperator4.Simplify seems only ensures the topology is correct
Integrate_management and RepairGeometry_management require feature classes
Eliminate requires Desktop Advanced which I don't have

My goal is to unionize the collection of IPolygon objects with no slivers?


texts: Interpolate vector layer and raster layer data into regular grid data using QGIS

texts:
I'm new to using GIS software. I obtained two different sets of data:

Topography data (raster, .asc)
Bathymetry data (shapefile, .shp, etc)

I include topography data using Add Raster Layer and bathymetry data using Add Vector Layer in QGIS. (See attached figure)

Is there any possible way to interpolate both different data sets like this? I have done some diggings, but they are mostly discussing how to interpolate vector layer data. I would like to get the water depth and land elevation in the same file after interpolation either in .xyz or .grd, so that I can use this as a computational domain for my simulation model.


texts:
I am trying to extract an ecoregion specified in the Ecoregions 2017 Resolve Data in Google Earth Engine in order to then to examine forest loss and gain (using the Hansen data set) from that ecoregion (tapajos xingu in Brazil). How do I specify the ecoregion tapajos xingu in the Ecoregions 2017 Resolve Data in Google Earth Engine?

var fc = ee.FeatureCollection('RESOLVE/ECOREGIONS/2017') .map(function(f) {
  var c = ee.String(f.get('COLOR')).slice(1);
  return f .set('R', ee.Number.parse(c.slice(0, 2), 16))
  .set('G', ee.Number.parse(c.slice(2, 4), 16)) .set('B', ee.Number.parse(c.slice(4, 6), 16));
});

var base = ee.Image(0).mask(0).toInt8();

Map.addLayer(base.paint(fc, 'R') .addBands(base.paint(fc, 'G') .addBands(base.paint(fc, 'B'))), {gamma: 0.3});

This is the code for Ecoregions 2017 Resolve Data


texts: Writing python script to batch add layers to multiple rasters?

texts:
So I have looked at existing tools and there is nothing to batch change symbology to a myriad of raster surfaces (20+). Basically my manual process is to right click on the rasters layer properties, symbology, show classified and import my previously defined color table layer onto the raster.

Is there a arcpy snippet I can use to do this batch wise to 20 plus rasters?

This is using ArcGIS 10.2.2


texts: adding photo to the pop up field in qgis2web?

texts:
I want to add charts on my web map in a way that when i click on the point the chart appears in the pop up field.

Do i have to add it as an image in the attribute?
Can i do that when the pictures are in my hard drive and not in a html form?


texts: Intersection of two sets in openstreetmap to find parent of any area

texts:
I am trying to determine what the parent administrative area of any given area. The idea is that the closest parent is the administrative area with the highest/deepest admin_level that is shared as a parent by all points of the area boundary. Her the pseudo-algorithm I would like to be able to reproduce with the OSM API.

For a given area:
For n in all nodes of its boundaries:

parents are the administrative area the node n lies in

potential_parent = intersect(list of parents)

unique_parent = parent in potential_parent where admin_level == min(admin_level)

Unfortunately I wasn't able to define the intersection of the different sets of parents for each node composing the boundary of the area of interest

Hehre is the solution I came up with:

rel(426391772) -> .obj;
.obj>;
is_in;
area._[boundary=administrative] -> .parents;
foreach .parents -> .shared(
  (.shared; - (.shared; - ._;););
);
out tags;

The for loop was my attempt to create the intersection set. But it doesn't return anything, while the .parents set is not empty.
I suspect the difference is not working as I thing it would.

The rest of the request to filter the good area is as follow and seems to work:

area.shared[boundary=administrative](if:t["admin_level"] == max(t["admin_level"]));

If you have a better approach to this problem, I am happy to hear it.

Here is an example of my request with the turbo-overpass.


texts: Locating date, scale, and coordinate system of shapefile in ArcMap?

texts:
I am trying to locate the date, scale, and coordinate system/projection of a shapefile in ArcMap but when I right click on the file in the table of contents and go to properties, I don't see the date or scale under the source tab.

I tried clicking on a feature with the identify tool but that also doesn't display the information.

How can I view this information from within ArcMap?


texts:
I try to replace part of text in 23 mxd files from "d" into 'map num: '. The code run without errors, but the result in the maps is only the character "m" instead of 'map num: '
I run this code:

import arcpy
from arcpy import env

env.workspace = r"G:\PROJECTS\road70-75"
counter = 0
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(
    r"G:\PROJECTS\road70-75\\" + mxdname)
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if elm.text.startswith(u'd'):
            elm.text = elm.text.replace(u'd', u'map num: ')[0]
            print 'elm.text','\n'
            counter = counter + 1
    mxd.save()
del mxd
print 'total maps: ', counter


texts: Visually correcting systematic shift between two geodatasets using ArcMap

texts:
I often have two datasets that are obviously not georeferenced correctly. Most of the time a simple translation (shift by x/y) could fix this.
Is there an easy way to do this e.g. by specifying one or two reference points in the reference raster and the raster that is to be corrected?
Is there a simple WYSIWYG-approach (e.g. in ArcMap) to manually/visually reference one dataset with some points to another and perform the transformation?
This seems like such a basic and straightforward task, that I think a tool should exist for it, but I know of none amongst the several tools at my disposal (ArcGIS, QGIS, R, Python, GDAL, ...).
I am aware of Shift Raster (Data Management) in ArcMap, but it requires known X/Y coordinates for the shift and doesn't work on vector-datasets. A workaround would be to extract coordinates of the two points, manually calculate the shift, and then use it with the tool. I'm hoping for something more elegant though.
I also found Transform Features (Editing) in ArcMap, which is fine but this too is restricted (only feature datasets) and circumstantial (as far as I understand one has to create features to represent the transformation).
E.g. for the data below I know that the green dataset is positioned right. Therefore I'd like to perform the translation represented by the yellow arrow (x -> o) on the complete grey raster data (and actually two other datasets as well). I'd love to be able to do this by just clicking at those two points and then on "shift".


texts: Renaming duplicates with any vacant value within a certain interval

texts:
The values of my column of interest (col1) should range from 0-999. I'm looking for a way to assign one of each duplicate a new value that lies within my this range 0-999.

Only features that are identical in both col1 and col2 are considered duplicates. Is there a way to identify "vacant" numbers and assign those to duplicated values.

At this point I use this code to identify duplicates and set new values for both of them:

UPDATE table AS t1
  SET col1 = t1.id || t1.col1
  WHERE EXISTS (
    SELECT *
    FROM table AS t2
    WHERE t1.id <> t2.id
      AND t1.col2 = t2.col2
      AND t1.col1 = t2.col1
  )

Any ideas?


texts: How to style Geofabrik shapefiles to resemble Open Street Map?

texts:
I have downloaded the shapefiles for the UK from http://www.geofabrik.de/data/download.html
I have loaded the shapefiles in to QGIS.
Could someone please explain to me or point me in the direction of how to correctly style the shapefiles to resemble Open Street Map, I have tried to google it and all the .qml files I find don't seem to do the job.


texts: Getting Google Maps "Areas of interest" Polygons

texts:
Is there a way to obtain the polygons corresponding to the the yellow-shaded "areas of interest" / areas of commercial activity from Google Maps?
Either through

a Google API
capturing the yellow zones in the tiles as polygons (raster to vector)
another data source (OpenStreetMap?) that also has "areas of interest" data
Something else


texts: How to preserve NoData areas when changing data type in GeoTiff with QGIS or gdal_translate?

texts:
I have a float geotiff raster with values 0, 1 and 2, along with a nodata area.  I'm trying to convert the file type do Byte but I'm not being able to preserve the nodata area. It gets converted to 0 and mixed with the data.

Here is the gdal_info of the source dataset

(output from gdalinfo clipped...)
Center      ( -50.5000000, -14.0000000) ( 50d30' 0.00"W, 14d 0' 0.00"S)
Band 1 Block=4701x1 Type=Float64, ColorInterp=Gray
  NoData Value=255

I then run:
gdal_translate -ot Byte teste.tif output.tif

The NoData value gets preserved, as gdalinfo returns

(output from gdalinfo cliped...)
Center      ( -50.5000000, -14.0000000) ( 50d30' 0.00"W, 14d 0' 0.00"S)
Band 1 Block=4701x1 Type=Byte, ColorInterp=Gray
  NoData Value=255

But looking at the image in QGis shows that the nodata area was transformed to 0

Here is the original raster in Float64 format

And here is what I get after gdal_translate

I'm using GDAL 2.2.2 in Ubuntu


texts: Using mapping/display tools in ModelBuilder?

texts: Does ModelBuilder have anything equivalent to the arcpy.mapping module, specifically zoomToSelectedFeatures, ExportToPDF, and TextElement?


texts:
I need to define a data structure which for convenience reasons for users using my software should be related to WGS84 - if possible. So that it makes it easy for users to use my software, convert from e.g. GPS data they get from Android or iOS system calls, and convert these forth and back to "my" coordinates.

However, my mapping application uses various UTM projections, depending which part of the world is shown. Thus I want to abstract away the need for a user to deal with various projections.

Ideally I just provide WGS84 as single projection - but I need height(!) since this is a 3D application and objects are on different floors with different elevation.

I've read some article here detailing about elevation here: https://geozoneblog.wordpress.com/2015/03/27/wgs84-versus-nad83-vertical/

What is the best way to provide WGS84 with elevation info, or shouldn't I do this?


texts: How to aggregate other columns when ST_ClusterIntersecting function is used

texts:
As stated in the PostGIS documentation

  ST_ClusterIntersecting is an aggregate function that returns an array
  of GeometryCollections, where each GeometryCollection represents an
  interconnected set of geometries.

Is it possible to aggregate other columns along with this geometry aggregation?
Particularly I would like to get arrays of IDs corresponding to the geometries in each aggregated geometry collection.


texts:
I am trying to generate a tile with Mapnik with a GeoJSON. The stylesheet.xml with the part to render the GeoJSON is the following:

<Style name="cities" filter-mode="first">
  <Rule>
    <MarkersSymbolizer width="50" fill="#ff4455" stroke="#881133" allow-overlap="true" ignore-placement="true" />
  </Rule>
</Style>
<Layer name="cities"
  srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs">
    <StyleName>cities</StyleName>
    <Datasource>
       <Parameter name="file"><![CDATA[/home/uli/projects/tiles/test.geojson]]></Parameter>
       <Parameter name="layer"><![CDATA[OGRGeoJSON]]></Parameter>
       <Parameter name="type"><![CDATA[ogr]]></Parameter>
    </Datasource>
  </Layer>

It does not render the GeoJSON. Also it doesn't throw an error. If I use a shapefile with the same settings (except type and file) it works fine. Any ideas why it is not rendered?

I used TileMill to generate the XML. The GeoJSON is shown properly in TileMill.

The GeoJSON file I used is https://gist.github.com/anonymous/2a140b557716c4237e96.


texts: Polygons not appearing when creating shapefile based on criteria using GDAL and Python (area calculated properly)

texts:
I am looping over the features of a shapefile containing multiple polygons. If  certain criteria are met, I want to save the current feature in a new shapefile.

I am doing this using GDAL 1.10.0 with Python 2.7.3

The code runs, creates the files I want, calculates the area correctly. In the shapefile that is created, the attribute table contains the propoer fields ('name' and FID). Two things are wrong:

No polygons are visible (in ArcMap). This seems to me particularly weird, because the polygons seem to me proper, in the sense that they are closed (and the area is calculated properly).
The coordinate system seems off (in Arc Catalog there appears in the bottom -20,9667 0,0000 (are these degrees??)), despite the fact that I am creating a projection file (they are in German national grid 3)

Any hints as to what the problems could be? I tried to put the relevant excerpt of the code below:

from osgeo import gdal, gdalnumeric, ogr, osr, gdalconst

### prepare shapefile into which to save
driver = ogr.GetDriverByName('ESRI Shapefile')
path_matched = r'.../.../path'
datasource = driver.CreateDataSource(path_matched)
lyr_matched = 'matched_%s' % cur_id  # from above
# Remove output shapefile if it already exists
name_shpfile = lyr_matched + '.shp'
if os.path.exists(os.path.join(path_matched, name_shpfile)):
    driver.DeleteDataSource(os.path.join(path_matched, name_shpfile))
layer = datasource.CreateLayer(lyr_matched, geom_type=ogr.wkbPolygon)
#layer.ResetReading()
field_name = ogr.FieldDefn("name", ogr.OFTString)
layer.CreateField(field_name)

### prepare file in which I want to write the areas of matched polygons
path_area = path_matched
file_area = 'matched_area_%s.txt' % (cur_id)
fobj_area = open(os.path.join(path_area, file_area), 'w')

### original shape file containing all polygons
shp = ogr.Open(shp_all)  # os fobj from above
lyr = shp.GetLayer()

# loop over polygons in original shameful
for FID in featList:
    feat = lyr.GetFeature(FID)

    # this is from the shapefile
    cur_name = feat.GetField(1)

    # ... work on conditions ...

    if condition = true:
        myPoly = ogr.Geometry(type=ogr.wkbPolygon)
        cur_poly = feat.GetGeometryRef()
        myPoly.AddGeometry(cur_poly)
        cur_spat_ref = cur_poly.GetSpatialReference()

        cur_area =  cur_poly.GetArea()
        cur_area_string = "%s \t %f \n" % (gw_name_EZG, cur_area)
        fobj_area.write(cur_area_string)

        featur_dfn = layer.GetLayerDefn()
        feature_abc = ogr.Feature(featur_dfn)
        feature_abc.SetGeometry(myPoly)
        feature_abc.SetField("name", cur_name)
        feature_abc.SetFID(FID)
        layer.CreateFeature(feature_abc)
        myPoly.Destroy()
        feature.Destroy()

# generate projection file
prj_filename = lyr_matched + '.prj'
file = open(os.path.join(path_matched, prj_filename), 'w')
targetSR = cur_spat_ref
targetSR.MorphToESRI()
file.write(targetSR.ExportToWkt())
file.close()

datasource.Destroy()


texts: Getting value from IRow is slower while edit mode is on?

texts:
When I access a value in an IRow using get_Value(index), it usually takes 0 to 1 ms. But when edit mode is on, the same operation takes 0 to 15 ms (most common values are 0 and 10 ms)! What could be the cause of this?

It doesn't matter if I use a filegeodatabase or an Oracle database (with SDE).

What we're using the get_Value for is our own ORM implementation (setting the Entity's property values to the ones within the IRow). Could there be a better way to do this?


texts: Updating field in feature class using ArcPy cursor?

texts:
I want to update a field in my feature class. I tried to use a da.SearchCursor (code below) but was quickly reminded that tuples are immutable... I ran the code but it failed with a TypeError: 'tuple' object does not support item assignment.

How do I most efficiently retrieve all rows from a feature class and update a field of each row and save it back to the feature class?

with arcpy.da.SearchCursor(fc, fldNames) as sc:
    for thing in sc:
        tmp = 0

        #
        # do stuff here
        #

        thing[idxField] = tmp


texts:
I am using AGS security and a 2.4 flex application.
I would like to login with a login screen like available in the samples.
Then utilize tokens for fc security on my services.

I see a sample for both of these on esri.resources.com but nothing for both together.
Has anyone done this that can share?


texts: Changing cardinality of existing relationship class in enterprise geodatabase?

texts:
Is it possible to change the cardinality of an existing relationship class in enterprise geodatabase?

I want to alter a 1-m relationship class to m-m


texts: Is there a well-known way to compare the relative accuracy of 2 latitude/longitude results?

texts:
I have 2 geocoding results (in latitue and longitude) for a given address from different services. I'm interested in describing the results in relation to one another in terms of relative accuracy and am wondering whether there's a well-known algorithm for that comparison in the GIS world.

Let's say we have Result A, which we assume for the sake of this effort is 100% accurate. Then we have Result B and want to know inaccurate it is in comparison to Result A. Result A has 15 digits of precision while Result B has only 6, but even at the 6-digit precision the results are different.

Does taking the difference between the two results and dividing it by the "100% accurate" result and regarding that figure as the relative accuracy in percent even make any sense in GIS data? If not is there something else that does?

For whatever it's worth, this is at least partially an exercise in curiosity, since I'm fairly sure that at the resolution we need both results will be close enough.


texts: Tool for making sure two vector layers are identical?

texts:
I have to constantly compare new project boundaries with old project boundaries (both are polygon shapefiles) to make sure that both the outside boundaries and the inside segments haven’t changed.

Is there a tool to check whether the two shapefiles are identical? So I don’t have to check each visually.


texts: Adding formatted sequence Number based on location?

texts:
I have Feature Point DATA. I need to add some sequence number to the attribute of each point, based on their location in the display map. In my attribute it'S ordered in the opposite way.

In attribute table its start from there..

I tried to descending the ObjectID first and put this code from ESRI but the result follows ObjectID sequence..

rec=0
def autoIncrement():
 global rec
 pStart = 1
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec += pInterval
 return rec

I tried to search how to do add sequence on descending way.
I found some advise from this but still cant do because it's too complicated for my data.
My data just contains ObjectID and shape.

So what I want to do is to create sequence number opposite way from ObjectID or create descending sequence number for the simple.

How do I do it?


texts:
I have thousands of pdf's that have meaningless names, and are scattered across several folders. I need to rename them, and fortunately most of them are referenced through links in spatial layers, where the attribute tables have meaningful information about the files and a column called 'hotlinks' has a directory path pointing to the pdf's.

Is there a method in QGIS that can follow the paths in the attribute tables, and rename the pdf's in that path based on the attributes in the layer? I will have a new field that contains the new name, eg. concat("Field1","Field2","Field3"), so I'd like to rename the pdf's based on that.

The catch...the pdf's are scattered around and are not in the same folder, so the script can't just point to one folder.


texts:
I'm facing problems with this plugin. I can not even add a toolbar. This is what I'm doing:

In index.html file I'm using this lines:

<script src="https://unpkg.com/leaflet@1.0.2/dist/leaflet.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/leaflet.draw/0.4.2/leaflet.draw.js"></script>

After, in the script file I'm using these lines to add the controls:

var drawnItems = new L.FeatureGroup(suc,crs);
map.addLayer(drawnItems);
var drawControl = new L.Control.Draw({
         edit: {
             featureGroup: drawnItems
         }
     });
     map.addControl(drawControl);

And I get this error:

TypeError: L.Control.Draw is not a constructor
http://xxxx.com/xxxxx/js/map.js Line 34

I've been searching in forums, the docs, and in several related pages but I haven't found any clue about what's happening.

If you need more information about it, just tell and I would be pleased to add it.


texts: Error displaying/cleaning a binary image in Google Earth Engine (from Sentinel-2 imagery)

texts:
I'm working on classifying the 'active' part of rivers using Sentinel-2 imagery in GEE (including water and the unvegetated, exposed sediment around the edges of the river).

The goal is to output a cleaned, binary image mask that represents the active river.

My code is working reasonably well and successfully displays the output layer over small areas (e.g. https://code.earthengine.google.com/3d709645a96ba017d157dacbb164f636), but fails to display the output layer when the geometry polygon increases beyond approx 10 km x 10 km (e.g. https://code.earthengine.google.com/25ce08a90fb05cec02cf488a8a455012).

The error message is:
"Output: Tile error: Output of image computation is too large (5 bands for 3260502 pixels = 124.4 MiB > 80.0 MiB)"

The error says that 5 bands are in the output image and wonder if this is partly causing the issue? However, printing the output binary image to the console shows that only a single band is present.

I'd hoped that the binary format would reduce the image size and allow for analysis over larger areas (without the need to break down into smaller tiles e.g 10 km by 10 km).

Any tips/workarounds would be appreciated.


texts: Calculating vegetation indices on pansharpened landsat 7 images

texts:
I have a problem with calculating vegetation indices like EVI, GSAVI, SAVI on Landsat 7 images after I done a Brovey Transform, Modified Resolution Merge and HCS transformation in Erdas 2013.

My study is based on comparing vegetation indices before and after pansharpening.

I didn't have any problems with calculating vegetation indices before pansharpening.

For calculating vegetation indices after pansharpening I first did TOA, and then with those images I performed pansharpening and then realised I need NIR band for calculating vegetation indices and in pansharpened images I only have R,G,B band.

What can I do?


texts: OpenLayers get closest feature over multiple vector layers

texts: My OpenLayers map has multiple, vector layers, each with unique features. I need a way to find the closest feature to a location across all the layers, in the same way that getClosestFeatureToCoordinate does for a single layer. I cannot find any OL method that would allow me to do determine that.


texts: Geoprocessor won't execute with ArcObjects

texts:
When I use the geoprocessor in my ArcObjects code for creating a feature in the sdeworkspace it doesn't work and returns an empty path for geoprocessor output.

How can I solve this problem?

        IWorkspace workspace = workspaceFactory.Open(connectionProps, 0);
        IFeatureWorkspace workspace  =(IFeatureWorkspace)workspaceFactory.Open(connectionProps, 0);
        IMxDocument mxd = ArcMap.Application.Document as IMxDocument;
        IMap map = mxd.FocusMap;
        ILayer layer= map.Layer[0] as ILayer;
        IFeatureLayer featurlayer= layer as IFeatureLayer;
        IFeatureClass fc = featurlayer.FeatureClass;
        IWorkspace wrk = workspace as IWorkspace;
        IDataset dd= outputFeatClass as IDataset;
        Geoprocessor gp = new Geoprocessor();
        FeatureClassToFeatureClass sss= new FeatureClassToFeatureClass();
        sss.in_features=fc;
        sss.out_path = wrk.PathName;
        sss.out_name=featurlayer.Name;
        gp.Execute(sss,null);


texts: Displaying percentages next to legend items in ArcMap?

texts:
I'm using ArcMap 10.4.1 and I have a map showing points on a map.

Each point represents a measurement of a specific radio frequency and has an attribute of signal strength.

The legend displays discrete bins of data within a range of possible radio frequency values.  I want to display the percentage of total points in each bin.

When I use the Legend Properties I'm able to display the count of points in each bin (checked 'show feature count' box in image below), but not the percentage.

This displays the count (in parenthesis) in each bin like in the map below:

How can I add percentage instead of count to my legend?


texts: Getting all the Raster pixel value from PostGIS

texts:
I just uploaded the raster using following code,

raster2pgsql -s 4326 -I -C -M raster.tif -F -t 100x100 public.raster_table | psql -U postgres -p 5432 -h localhost -d raster

It takes around 15 min to upload to database. After complete the process, I checked the table using SELECT * FROM public.raster_table;. It has lots of  rows (around 4000) with serial rid. I want to view the each pixel value using ST_value. I already checked the documentation. In documentation, It provides following method to show all the pixel value.

SELECT x, y, ST_Value(rast, 1, x, y) As raster_val,
FROM raster_table CROSS JOIN
generate_series(1, 1000) As x CROSS JOIN generate_series(1, 1000) As y
WHERE rid =  2 AND x <= ST_Width(rast) AND y <= ST_Height(rast);

But I am confused with WHERE condition. In example it only use rid=2. But I have lots of rid. Is this condition really checks the all the pixel value?


texts:
One of the databases I use is throwing a DBMS error when editors are editing the database through ArcMap. It says that the layer was unloaded due to memory pressure. This error doesn't occur for every editor when they are working with the versioned layers in an editing session. It seems to be an intermittent issue limited to one database when editors are adding data to feature classes.
I have no idea what that means. Can anyone help?


texts: How to indicate a range of zoom on the map?

texts:
Basically what I need and I want to know is how you can do to tell the map that respects me a zoom range to display the data, ie, I want to see the map from afar without information, and as the zooming you can go see the information.

Can anyone help me with this?


texts: Dividing region in 9 sub-regions in Google Earth Engine

texts:
I am studying the spectral response of PV plants in GEE. I need to divide a polygon (geometry), which is shaped like the plant, in 9 sub-polygons, to investigate how each corner of the plant reacts.
Is there any way I could not do it by hand, and let GEE do it?


texts: Correcting Worldview-3 data to ground reflection/bottom of atmosphere with open source?

texts:
I have some worldview-3 imagery that needs to be corrected from digital numbers to actual reflection values (preferably bottom-of-atmosphere). But the tricky part is I am looking for open source options.
I see that this is possible for worldview-2 with the Orfeo toolbox but worldview 3 isn't supported. Is there a way to work around this by chance or does anyone have any other suggestions on python libraries/QGIS tools, etc. ?


texts: arcpy Geometry tools give incorrect results?

texts:
I'm new to python and arcpy but I've been trying to learn how to run geoprocessing tools using arcpy. I have (for this example) 2 shapefiles, a and b, which are both polygons. From plotting them in arcMap, I know that b is found entirely within a.

I would like to write a script that will take these two shapefiles and give me a True/False response about whether or not one polygon contains the other. I was thinking of using the arcpy.Geometry.contains() function.

I have the following:

>>> a_shp="C:\\Data\\shapefile_a.shp"
>>> b_shp="C:\\Data\\shapefile_b.shp"

>>> #read in as geometry
>>> a=arcpy.CopyFeatures_management(a_shp, arcpy.Geometry())
>>> b=arcpy.CopyFeatures_management(b_shp, arcpy.Geometry())

>>> arcpy.Geometry.contains(a[0],b[0])
False
>>> arcpy.Geometry.contains(b[0],a[0]) #maybe I got the order wrong? Trying the other order
False

>>> #Maybe I'm getting the syntax wrong, trying something else
>>> a[0].contains(b[0])
False
>>> b[0].contains(a[0])
False

I get no error messages but clearly something isn't working properly. Am I doing this right?


texts: How to restrict user access to map services on ArcGIS for Server?

texts: If we have all our map services published on one server, can we configure it such that our users can only see or access a only a subset of the map services?


texts: An unhandled exception of type System.AccessViolationException

texts:
I am trying to upgrade a desktop extension from ArcGIS 10.0 to ArcGIS 10.2. I am getting this error:

An unhandled exception of type 'System.AccessViolationException' occurred in GIS.Metadata.dll.

Not sure what may be causing this. My Environment: VS2013, ArcGIS 10.2, VB.NET, Windows 7 (64 bit). I have also added code where it’s failing.

pGxCatalog = New ESRI.ArcGIS.Catalog.GxCatalog
'// connect to root ...
pGxCatalog.ConnectFolder(IO.Path.GetPathRoot(pPath))

'// connect to folder via Location property
pGxCatalog.Location = IO.Path.GetDirectoryName(pPath) '

Fails here

pCatObject = pGxCatalog.GetObjectFromFullName(pPath, numFound)

IF I set the location manually (say C:) it errors on the above line.


texts:
I have 190 layers in shape file with attribute filled in lower case. I want to change data of all tables in uppercase using Qgis.

I am able do to do it one by one using upper in expression.

How to do for all layers in one go?


texts: Renaming bands using Google Earth Engine JavaScript API

texts:
I would like to rename the band "B4_stdDev" to "Red_stdDev", and the band "B10_stdDev" to "Temperature_stdDev".
Example code:
var b4_stdDev = ee.Image().rename('B4_stdDev');
var b10_stdDev = ee.Image().rename('B10_stdDev');

var img = ee.Image.cat(b4_stdDev, b10_stdDev);
print(img);

This is a simple example. I actually have an image with approximately 60 bands with names like "B2_mean", "B2_savg", "B2_ent", etc.


texts: Geocoding several addresses to Multipoint geometry

texts:
I have a csv-file of different companies, see sketch below. Each company has at least one address. Some can have more than 50 addresses. There is one address per cell. In total my csv-file has around 3000 rows.

+--------------+------------------+------------------+------------------+
| company_name |    address_1     |    address_2     |    address_3     |
+--------------+------------------+------------------+------------------+
| name_1       | some_address_1   | some_address_2   |                  |
| name_2       | 2_some_address_1 |                  |                  |
| name_3       | 3_some_address_1 | 3_some_address_2 | 3_some_address_3 |
+--------------+------------------+------------------+------------------+

My goal is to geocode these addresses to points and convert them to a Multipoint geometry. And then represent each company by a single entity with its different locations.

I have tried with MMQGIS and Google Earth but I don't know how to convert other addresses.

My main trouble is about geocoding several columns in a same row and merge the points to a multipoint geometry.

And I would like to add a new column containing a multipoint geometry with the geocoding of the cells which are not empty.

Is there a solution with geopandas or PostGIS?

The thing is that once the process is completed, I intend to store the result in a PostGIS database and display it on QGIS.


texts: Getting ERROR 999999 when trying to use arcpy.UploadServiceDefinition_server()?

texts:
Here is my code

def upload(name, client, working_dir):
  arcpy.env.overwriteOutput = True
  arcpy.SignInToPortal('https://www.arcgis.com', username, password)
  outdir = working_dir + "\\PROJECTS\\" + name + "\\sd\\"
  aprx = arcpy.mp.ArcGISProject(working_dir + "\\PROJECTS\\data_holder\\data_holder.aprx")
  m = aprx.listMaps()[0]
  print("Uploading data for " + name)
  lyrs = m.listLayers()
  arcpy.mp.CreateWebLayerSDDraft(lyrs, outdir + name + ".sddraft", name, "MY_HOSTED_SERVICES", "FEATURE_ACCESS", name, True, True, False, False, False, name, name, name, None, None )
  arcpy.StageService_server(outdir + name + ".sddraft", outdir + name + ".sd", 102)
  arcpy.UploadServiceDefinition_server(outdir + name + ".sd", 'My Hosted Services')

This is the error I am getting. I have tried multiple different ways to put in parameters.

  Uploading data for test_project

  Traceback (most recent call last):

  File "main.py", line 304, in module

  upload(name, client, working_dir)

  File "main.py", line 242, in upload

  arcpy.UploadServiceDefinition_server(outdir + name + ".sd", 'My Hosted Services')

  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\server.py", line 954, in UploadServiceDefinition
      raise e

  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\server.py", line 951, in UploadServiceDefinition
      retval = convertArcObjectToPythonObject(gp.UploadServiceDefinition_server(*gp_fixargs((in_sd_file, in_server, in_service_name, in_cluster, in_folder_type, in_folder, in_startupType, in_override, in_my_contents, in_public, in_organization, in_groups), True)))

  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing_base.py", line 506, in
      return lambda *args: val(*gp_fixargs(args, True))

  arcgisscripting.ExecuteError: ERROR 999999: Error executing function.

  ERROR: code:400, Invalid URL, Invalid URL, Bad syntax in request.

  Failed to execute (UploadServiceDefinition).


texts: How to identify concentrated areas with multiple point values to produce a "heat map" appearance

texts:
I have roughly 1.6 million points statewide with values ranging from 1 - 4. I'm looking to identify concentrations of points based on their values (1,2,3,4). I know there are several different routes you could take on this approach, but I am struggling to find one that produces a visually clean output. I am wanting to produce a map that would look similar to a "heat map." Here are a few of the different methods I've tried so far.

1) Converting all the points into a raster and using the most frequent cell assignment. The output is what I would expect, but I'm not a fan of the pixelated look. I used the smoothing tool and that cleaned up the raster cell edges, but I'm still not truly satisfied with the look of the end product.

2)The second approach I've tried is neighborhood point statistics using majority as the statistical type. This approach seems very subjective depending on the cell size and neighborhood size. The output however appears to be more of what I was hoping the end product would look like. However when looking at the individual point data, some of the output looks questionable.

Can anyone think of a more appropriate approach that would allow for more statistical analysis with the end product looking similar to a heat map? I've been looking into the hot spot analysis tool. I don't think that tool will work with my data since I have multiple point values. I am using ArcGIS 10.2 info.


texts: Google Earth Engine: get data ready for ML export

texts:
I need to export data from GEE to get processed in MATLAB. I want to export CSV files of a region, aggregating meteo data ( CAMS, ERA5, MODIS AOT) with sentinel-2 images. These images have already been processed to have 12 images, where each one is the mean over a month. I wanna reduce this images to a single number,with a mean over a region (called plant), then construct an array, then export it to a csv file.
Each row of the file is a month, each column is a feature (meteo data) or an output (sentinel-2 bands). Then I will import this data for further processing in MATLAB.
I know that I have to use reducers for this matter, like ee.Reducer.mean(),but how exacly can I obtain such a structure, and what is the right reducer to get a spatial mean, and not a temporal one?
So far, this is my code: https://code.earthengine.google.com/874acca5c75fd2cc6135c3c7055ba8e9


texts: Large scale map projection (ca. 20 km) in pyproj results in wrong distance

texts:
I would like to transform from lon/lat (EPSG:4326) to a local x/y coordinate system for large scale equidistant maps with a map size of ca. 20 km. I need to do this mapping all over the globe apart from regions close to the poles.
Looking at https://projectionwizard.org/ it suggests for 1x1 degs something like +proj=cass +lon_0=-121.5 +datum=WGS84 +units=m +no_defs.
However, when looking at distances e.g. high up in Norway, the distances are way off. I must be doing something wrong, but I have no clue what it is.
self.trans_proj = pyproj.Transformer.from_crs(
  'EPSG:4326',
  {
    'proj': 'cass',  # tmerc
    'x_0': '500000.0',
    'lon_0': str(self._lon),
    'lat_0': str(self._lat),
    'k_0': '0.9999',
    'datum': 'WGS84'
  }
)

x, y, = self.trans_proj.transform(lon, lat, radians=False)

I set the lon_0 and lat_0 attributes to the current lon/lat of the centre of the map.


texts:
Is it possible to lock a QGIS project file, so only I can open it with a password?

I work with a QGIS 2.18 version and I'm looking for a solution inside the QGIS software. in my map, i use only with regular layers and use not PostGIS.
I red Password Protect QGIS project and it didn't answered my question.


texts: Joining wind speed and wind direction using Google Earth Engine

texts:
I have successfully exported wind speed to a csv file, and am trying to add wind direction. I know my wind direction calculation works and I can see it when printing ERA5total. If I understand correctly, I need to use an inner join, but I only end up with an empty csv file.
//  1. Setup
var countries = ee.FeatureCollection("FAO/GAUL/2015/level0");
var NL = countries.filterMetadata('ADM0_NAME', 'equals', 'Netherlands');
var start_period = ee.Date('2018-07-01');
var end_period = ee.Date(new Date().getTime());

var ERA5 = ee.ImageCollection("ECMWF/ERA5_LAND/HOURLY")
     .filter(ee.Filter.date(start_period, end_period))
     .filterMetadata('hour', 'equals', 14);

var ERA5NL=ERA5.map(function(image){
   return image.clip(NL)});

//  2. Calculate wind speed as windspeed band
var ERA5windspeed = ERA5NL.map(function(image){
  var wind_10m = image.expression(
    'sqrt(u**2 + v**2)', {
      'u': image.select('u_component_of_wind_10m'),
      'v': image.select('v_component_of_wind_10m')
  }).rename('windspeed');
  var time = image.get('system:time_start');
  return wind_10m.set('system:time_start', time) } );

//  3. Calculate mean wind speed for entire NL per day
var ERA5meanspeed = ERA5windspeed.map(function(image){
  var meandict = image.select('windspeed').reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: NL
  })
  return image.set(meandict);
});

  //  4. Calculate wind direction as winddir band
var ERA5winddir = ERA5NL.map(function(image){
  var winddir = image.expression(
    'mod(180 + (180/3.14) * atan2(v,u),360)', {
      'u': image.select('u_component_of_wind_10m'),
      'v': image.select('v_component_of_wind_10m')
  }).rename('winddir');
  var time = image.get('system:time_start');
  return winddir.set('system:time_start', time) } );

//  5. Calculate mean wind direction for entire NL per day
var ERA5meandir = ERA5winddir.map(function(image){
  var meandict = image.select('winddir').reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: NL
  })
  return image.set(meandict);
});

//6. Join wind speed and wind direction datasets

// Use an equals filter to define how the collections match.
var filter = ee.Filter.equals({
  leftField: 'system:index',
  rightField: 'system:index',
});

 // Define the join.
var innerJoin = ee.Join.inner('primary', 'secondary');

// Apply the join.
var ERA5total = innerJoin.apply(ERA5meanspeed, ERA5meandir, filter);

// Display the result.
print('Join', ERA5total);

//7. Filter out wind speeds > 3 m/s and keep only yyyymmdd in output table
var ERA5x = ERA5total.filter(ee.Filter.lt('windspeed',3.0))  //select only wind speeds < 3 m/s
  .map(function(image){
      var yyyymmdd = ee.String(image.get('system:index')).slice(0,8) //slice date: only yyyymmdd left
    return image.set('date', yyyymmdd)
});

Export.table.toDrive({
  collection: ERA5x.sort('windspeed'),
  description: 'NL_mean_wind',
  selectors: ['date', 'windspeed', 'winddir']
});


texts: Selecting features within polygon from another layer using R or Python?

texts:
For Selecting features within polygon from another layer using QGIS? the solution was to:

  use the "Vector->Research tools->Select by Location"

However, I would like to know how to do the same thing with either an R or Python script.


texts: How to delete points with the same information (housenumber) but different position within polygons (buildings) and keeping just one in QGIS?

texts: I´m working with shape files. I have one point layer with house numbers, which I want to de shown over the polygon layer buildings. My Problem is that almost every builing has got 2 or more points with the same number. I want to keep just one. But there are also points within some buildings with different information e.g. house number 2a, 2b, 2c and this should stay like it is.


texts: Random selection of points in SpatialPointsDataFrame R object with distance constraint

texts:
I need to select points randomly (without replacement) from a SpatialPointsDataFrame (input object) in R with a minimum distance of 3,000 meters between them. I want to get random 50% of all points. I know that "Create Random Points" tool from ArcGIS, as mentioned before Randomly sampling points in R with minimum distance constraint, can do this processing, but I really need
to do this inside R. I tried to use sample() function but I still did not realised how to set the geographical constraint. I tried to run QGIS inside R, but it seems that QGIS does not have a tool for that.

> input
class       : SpatialPointsDataFrame
features    : 205
extent      : 203294.5, 259880.6, 7600123, 7668676  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=23 +south +ellps=GRS80 +units=m +no_defs
variables   : 3
names       :     Sites, Long_X,   Lat_Y
min values  : CPF - 050, 203295, 7600120
max values  : JES - S94, 259881, 7668680

> head(input)
Sites Long_X   Lat_Y
1 CPF - 050 235441 7617150
2 CPF - 052 234106 7615740
3 CPF - 054 232683 7614280
4 CPF - 056 233863 7614420
5 CPF - 058 234012 7612890
6 CPF - 062 236929 7612850

rdm_input <- sample(x= nrow(input), size=(0.5*205), replace= FALSE)
[1] 167 189 79 80 126 129 144 100 4 109 170 72 123 73 132 93 169 5 134 176 196
158 152 183 23 136 180
[28] 130 12 142 179 11 13 66 22 2 96 29 54 137 120 171 184 36 113 3 81 115  30
85 61 162 98 102
[55] 103 181 90 133 56 174 76 201 150 197 14 86 118 121 28 97 160 178 1 186
141 163 172 32 65 168 74
[82] 114 182 128 131 67 70 165 187 185 69 17 194 16 154 119 192 156 106  25
101 198 105 108 151 125 190 10
[109] 84 38 51 161 40 94 45 19 145 75 42 122 117 49 87


texts:
I'm trying to automate adding styles to layers, with each layer having multiple styles using pyqgis.

Couldn't find any answers except this unanswered question:

Adding multiple styles to single layer with pyQGIS?

I did manage to add multiple styles using QgsMapLayerStyleManager:

layer = iface.activeLayer()
layer.name()
mlsm = QgsMapLayerStyleManager(layer)
layer.loadNamedStyle('~/Styles/BasinAnalysis/Area/Total.qml')
mlsm.renameStyle("default","Area-Total")
mlsm.addStyleFromLayer("Area-Freshwater")
layer.loadNamedStyle('~/Styles/BasinAnalysis/Area/Freshwater.qml')
mlsm.addStyleFromLayer("Area-Unirrigated")
layer.loadNamedStyle('~/Styles/BasinAnalysis/Area/Unirrigated.qml')
etc...

It seems like the multiple styles are associated to layer, and I can also see the changes applies to layer:

mlsm.styles()
>>> ['Area-Freshwater', 'Area-Other', 'Area-Saline', 'Area-TWW', 'Area-Total', 'Area-Unirrigated']

The problem is that styles don't appear when I right click on Layer -> Styles, it only shows the default layer.

I'm guessing I need to associate the styles from the StyleManager back to the layer, but not sure how.


texts: Adding z values to contours without attribute data

texts:
I have a set of contour lines with 5 meter interval which doesn't have any height set in the attribute data.

Now I'd like to assign (preferably automatic) z values to these contours since I want to change the interval to 10 meters. Does anyone know how to do this using open source tools?

Almost ten years ago I found an application where I could assign height to contours by dragging a line through some contour lines and tell the app what the interval was and the height of the first contour line. I don't remember what the app was called but that one would help me even though it would be nicer to add elevation info from a dem or something.


texts: Converting CSV file to GeoJSON while preserving data types?

texts: Assuming said csv file is properly structured with Latitude and Longitude.  I'm having problems find converters that won't mess up the data types (putting string around integers).


texts: Get GPS Coordinates on mouse tap/click on ArcGIS Esri Map in Windows Runtime Apps (Windows 8.1)

texts:
How to get the GPS coordinates (longitude, latitude) of the clicked/tapped point in ArcGIS Esri map?

I tried MapViewTapped event like:

private void myMapView_MapViewTapped(object sender, MapViewInputEventAtgs e)
{
  var x = e.Position.X; //e.Location.X;
  var y = e.Position.Y; //e.Location.Y;
}

Both gives large unexpected values.


texts: Downloading data directly onto Google Earth Engine (GEE)

texts: Is there a way to download external files (.tiff or .csv for example) directly onto my earth engine account?  Currently, to import external files, I have to download it first onto  my own computer, then upload it to the GEE.


texts: Error during load osm pbf file from QGIS

texts:
I was trying to print out Stamen map (Quickmap plugin) in very high quality. However, I was unable to do so, printed out results looks like chequer board or empty.

So, I was following the information from
How to create high quality map with QGIS and Stamen tiles and
https://anitagraser.com/2014/05/31/a-guide-to-googlemaps-like-maps-with-osm-in-qgis/

And I downloaded osm.pbf for asia from http://download.geofabrik.de/index.html here.

Problem is, OSGeo4W shell makes warning messages while I tried to transform pbf to db file:

C:\>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES D:\QGIS_Korea\test2.db D:\QGIS_Korea\asia-latest.osm.pbf
0...10...20...30...40...50...60...70...80...90...Warning 1: organizePolygons() received an unexpected geometry.  Either a polygon with interior rings, or a polygon with less than 4 points, or a non-Polygon geometry.  Return arguments as a collection.
Warning 1: organizePolygons() received an unexpected geometry.  Either a polygon with interior rings, or a polygon with less than 4 points, or a non-Polygon geometry.  Return arguments as a collection.
Warning 1: organizePolygons() received an unexpected geometry.  Either a polygon with interior rings, or a polygon with less than 4 points, or a non-Polygon geometry.  Return arguments as a collection.
100 - done.
Warning 1: organizePolygons() received an unexpected geometry.  Either a polygon with interior rings, or a polygon with less than 4 points, or a non-Polygon geometry.  Return arguments as a collection.
Warning 1: organizePolygons() received an unexpected geometry.  Either a polygon with interior rings, or a polygon with less than 4 points, or a non-Polygon geometry.  Return arguments as a collection.

Anyway it finished the conversion. But when I try to add them to my project via add layer -> Spatialite layer -> new -> choose file -> connect, it prints out following error message.

Spatial metadata check failed
Failure getting table metadata... is ****.db file really a SpatiaLite database?
table info on geometry_columns failed
database is locked

How can I escape from this situation?


texts: Adding Logo and Title to an HTML webmap created by Qgis2Web

texts: I would like to add a company logo stored in my machine to a webmap created with qgis2web in QGIS. Also, it would be nice if I can add a title as well. Anyone could help me with that?


texts: How to debug 'Layer is not valid error' from Python plugin?

texts:
I'm trying to understand what goes wrong with a QGIS python plugin.
Basically I'm trying to render a series of POINT.

The query works fine, at least giving it by hand from inside PostGIS DB:

"SELECT  gid, CAST (gid as text)|| '_' || (ST_DumpPoints(the_geom)).path[2] as key,  (ST_DumpPoints(the_geom)).geom as points FROM ways"

I concat gid and path[2] because after splitting the_geom (MULTILINE) field into several POINT, gid are duplicated and I think they are no longer usable as key (with this trick key should be unique, but it seems not solve my problem).

I tried to create the layer this way:

def createPointLayer(self):
    uri = self.db.getURI()
    query = "SELECT  gid, CAST (gid as text)|| '_' || (ST_DumpPoints(the_geom)).path[2] as key,  (ST_DumpPoints(the_geom)).geom as points FROM ways"

    uri.setDataSource("", "(" + query + ")", "points","","key")
    layerName = "NodeLayer"
    aLayer = self.iface.addVectorLayer(uri.uri(), layerName, self.db.getProviderName())

I think the query works fine (QGIS takes a long time before display the error. I guess something goes wrong after the query has been executed).

When the plugins run QGIS output the following error:

  Layer is not valid:

  The layer dbname='pgrouting' user='user' password='user' key='gid'
  table="(SELECT  gid, CAST (gid as text)|| '_' ||
  (ST_DumpPoints(the_geom)).path[2] as key,
  (ST_DumpPoints(the_geom)).geom as points FROM ways)" (points) sql= is
  not a valid layer and can not be added to the map

I already had to struggle against the Layer is not valid error, but I don't know how to debug it.
So I have 2 questions:

Could anyone tell me what's my mistake in this specific case?
How can I debug 'Layer not valid error'? How can I get more information on which is the problem?

EDIT:

QGIS log console display the following message:

  No key field for query given


texts: How to replicate ST_Length using python without saving to DB?

texts:
I'l like to be able to calculate the length in meters of a LineString, in python. I came across this question: Get LineString length in meters (Python, GEODjango)
but the recommended method seems to disagree with the calculation given by ST_Length (spatialite in this example) - which also agrees the actual length of the path.

I don't want to use ST_Length, I need to find the length in python. But I can't figure out why these calculations differ. To exemplify:

My Model:

class TestRoute(models.Model):
    path = models.LineStringField(null=True)

Using GeoDjango's Length annotation after creating an instance with an example path returns the correct length in meters:

In [1]: TestRoute.objects.annotate(length=Length('path'))[0].length.m
Out[1]: 41345.49910715875

But using the same instance:

In [2]: path = TestRoute.objects.all()[0].path
In [3]: path.srid
Out[3]: 4326
In [4]: path.length
Out[4]: 0.416972813296872
In [5]: path.transform(3857)
In [6]: path.length
Out[6]: 51830.398931961485

The difference is about 10.5k meters.

How can I get the total length of a LineString without needing to save it to a model?


texts: Adding point to existing attribute record, rather than editing attributes of existing point?

texts:
I'm using ArcGIS 10 to create point data. Normally I would start an edit session, select point as a construction tool, click where I want a point to be, and then go in and edit the attributes.

Say I forgot to do that and just went into the attribute table and added a few records.

Am I stuck deleting the records I've just added and starting over, creating the points first this time, or is there a way to associate these records with new points?


texts: Is there any online tool for converting coordinates between different EPSG?

texts: Every time I want to do conversion, I need to get into the console to run proj.4. Is there any online tools for converting coordinates between different EPSG?


texts:
I'm new to Python and trying to automate copying specified features from a bunch of nested geodatabases into a new geodatabase- script seems to work ok, except for the copying features part (the 2 lines commented out below) .

import os
import arcpy

top_folder = r"W:\Gisdata\Fish\Proj\Bathymetry"
output = r"D:\Fish\PROJ\Bathy\test\testmrg.gdb"
bathshore = []

for path, dirs, files in os.walk(top_folder):
    for d in dirs:
        if not d.endswith(".gdb"):
            continue
        gdb_path = os.path.join(path, d)

        arcpy.env.workspace = gdb_path
        all_fcs = arcpy.ListFeatureClasses()
        for fds in arcpy.ListDatasets('','feature'):
            for fc in arcpy.ListFeatureClasses('','',fds):
                all_fcs.append(fc)

        for a in all_fcs:
            if '_shore' in a:
##                arcpy.CopyFeatures_management(a, output+"sh")
##                bathshore.append(a)
                print "Exported %s" % a

print "Exported the following feature classes %s" % bathshore


texts:
I am a total PostGIS noob, and am having trouble understanding the finer points of SQL syntax. I am trying to take a large set of line data and clip it with a polygon. I then want to create a new table with the lines that are in the polygon I used to clip. When I run this I receive ERROR: syntax error at or near "ST_Split" I am not sure what is wrong with my syntax. I was hopping that some could explain to me what is wrong with this line of code. Additionally I hoping that some one could explain the looping and variable syntax. Coming form a python background it seems very complicated.

CREATE TABLE roads_clipped AS

ST_Split(
    (SELECT geom
    FROM roads),
    (SELECT st_union
    FROM iberian_peninsula)
);

CREATE TABLE roads_of_the_iberian_peninsula;

BEGIN
    FOR row IN
        SELECT geom FROM roads_clipped
    LOOP
        IF ST_Contains(
                (SELECT st_union
                FROM iberian_peninsula)
                ,row)
            THEN INSERT INTO roads_of_the_iberian_peninsula
                VALUES row;
        END IF;
    END LOOP;
END;


texts: how to Input polygon featuer using Rubber poly (ArcObjects)

texts:
Does any one know how to input polygon feature using Rubber poly?
I would like to know how to input Polygon featurer using rubber poly..
as of now i can input polygon feature. but some times , polygon are is negative value and can not show the label etc..

I assume i need to something about geometry setting.. but i do not know how to do it...

protected override void OnMouseDown(ESRI.ArcGIS.Desktop.AddIns.Tool.MouseEventArgs arg)
{
    base.OnMouseDown(arg);
    IMxDocument pMxDocument = (IMxDocument)ArcMap.Application.Document;
    IMap pMap = pMxDocument.FocusMap;

    IRubberBand pRubber = new RubberPolygon();
    IPolygon pPolygon = pRubber.TrackNew(pMxDocument.ActiveView.ScreenDisplay, null) as IPolygon;

    IFeatureLayer pFLayer = pMap.Layer[0] as IFeatureLayer;
    IWorkspaceEdit pWorkspaceEdit = ((IDataset)(pFLayer.FeatureClass)).Workspace as IWorkspaceEdit;

    pWorkspaceEdit.StartEditing(false);
    pWorkspaceEdit.StartEditOperation();

    IFeature pFeature = pFLayer.FeatureClass.CreateFeature();
    pFeature.Shape = pPolygon;

    pFeature.Store();

    pWorkspaceEdit.StopEditOperation();
    pWorkspaceEdit.StopEditing(true);
    pMxDocument.ActiveView.Refresh();
}


texts: Fuzzy least-coarse common partition algorithm

texts:
Given two different partitions of a shape (for argument's sake, two different administrative divisions of a country), how can I find a new partition into which both of those partitions fit, allowing for (and optimising) some error?

For example, ignoring the error, I want an algorithm that does this:

Perhaps it helps to express this in set terms. Using the following numbering:

I can express the partitions above as:

A = {{1},{2},{3,4,7,8},{5},{6},{9,10,13,14},{11},{12},{15},{16}}

B = {{1,2,5,6},{3},{4},{7},{8},{9},{10},{13},{14},{11,15},{12,16}}

A dot B = {{1,2,5,6},{3,4,7,8},{9,10,13,14},{11,15},{12,16}}

and the algorithm for producing A dot B seems straightforward (something like, if two elements are in a partition together in A (B) merge the partitions they are in in B (A) - repeat until A and B are equal).

But now imagine that some of these lines are slightly different between the two partitions, so that this perfect answer is not possible, and instead I want the optimal answer subject to minimizing some error criterion.

Take a new example:

Here in the left column we have two partitions with no common lines (apart from the outside border itself). The only possible solution of the above kind is the trivial one, the right column. But if we allow "fuzzy" solutions, then the middle column may be permissible, with say 5% of the total area being contested (i.e. allocated to a different subarea in each coarsened partition). So we might describe the middle column as representing the "least-coarse common partition with <= 5% error".

Whether the actual answer is then the partition in top row, middle column or middle row, middle column - or something in between, is less important.


texts: Selecting Different points in 2 different layers

texts:
I'd like to make an analysis on points.

I've 2 layers they have house numbers as an attribute. What i'd like to do is compare these 2 layers and select the different ones.

Geometrically there are average 3 meters difference between 2 layers. First i used select by location on arcmap and i tried to select the points which are in 5 meters radius then inverted the selection.But results are not satisfying. now i am searching for another method.

Do you have any ideas? I don't have any expierence in python or something else. I want to do that analyis with ArcMap tools.


texts: Displaying map in p.mapper framework for UMN MapServer?

texts:
I want to display map in p.mapper. I've done some goggling and looked at some examples but can't seem to make what I've seen work. Can someone recommend a link or example I can look at?

I'm using mapserver and p.mapper


texts: Merging thousands of rasters without downloading images

texts:
My AOI has about 8262 1m DEM files (7 GB size). I want to merge these files without downloading these images from https://elevation.fsdf.org.au/ using gdal.
Is there a way to merge thousands of rasters without downloading the images?


texts: How to deal with QGIS Time Manager error: no attribute 'animationFrameCounter'?

texts:
I am using QGIS 2.6 and using the time manager I get the error shown in the attached image.

Does anyone know what the problem is and could offer a solution?

!


texts: Trilateration algorithm for n amount of points

texts:
I need to find algorithm that can calculate centroid A (aka gravity center, geometric center, center of mass) from the figure where circles T1,T2,T3,T4,T5,..,Tn intersect
AND length of line R from centroid to farthest corner of mentioned figure

Following information is given:

T1   Latitude    =   56.999883   Longitude   =   24.144473   Radius  =   943
T2   Latitude    =   57.005352   Longitude   =   24.151168   Radius  =   857
T3   Latitude    =   57.005352   Longitude   =   24.163356   Radius  =   714
T4   Latitude    =   56.999042   Longitude   =   24.168506   Radius  =   714
T5   Latitude    =   56.994226   Longitude   =   24.15709    Radius  =   771

Result should look like this:
 A Latitude = XX.XXXXXXX    Longitude = XX.XXXXXXX  Radius = XX

As you probably already figured out, I am working on software that can find device location by closest Wifi Access Points or Mobile Base stations, as number of access points or base stations might change, I need an algorithm that can adapt to uncertain amount of points.

There are some similar questions here and here, but none of them exactly answers to my question.


texts: Filling holes between multiple polygons in ArcGIS Desktop

texts:
How can I fill holes or gaps between multiple polygons in ArcGIS automatically?
For example, creating polygon number 2 below.


texts: QGIS/Spatialite cluster analysis with density

texts:
beginner here
I am trying to analyse a large set of points (2 millions) over a set of polygons (3000), and particularly find clusters of points. Doing it in QGIS was taking too long/crashing the compputer so I uploaded my layers on a Spatialite DB and use the interface provided by QGIS.

I have managed to do what I want using a few scripts in a row, but the cluster analysis is based on measuring the area of the dissolved cluster of buffered points, which doesn't take into account density. I have tried to add "weight" to the points and do a count but I can't get the query to work because the query results in a multi-part polygon.

Any idea to simplify or improve this process to a lower number of steps?

Query to calculate the number of points within the polygons:

SELECT b.name AS NameID,
count(a.geom) AS density
FROM
Point_layer AS a,
Polygon_layer AS b
WHERE ST_Intersects(a.geom, b.geom) = 1
group by NameID

Query to get the ID of points that intersects polygons (used to reduce the size of the point database to query:

 SELECT
b.Name AS Name_id,
a.PointID AS Point_ID
FROM
Point_layer AS a,
Polygon_layer AS b
WHERE ST_Intersects(a.geom, b.geom) = 1
group by Point_ID

Then I sort out the points I need to a new layer, add a buffer and dissolve:

Select a.id, (st_unaryunion(st_collect(ST_Buffer(a.geom, 7)))) as geometry
FROM point_layer_filtered as a

Then I split the resulting polygon into multipart into singlepart, and calculate area.

What I would like is to use something else than ST_unaryunion that allow to retain the number of points per clusters..

Is there a more efficient strategy that doesn't require buffers? Something that just group points together?
I guess I could rerun the first script on the dissolved buffers but that would take forever

EDIT, Solution:

Response that might help others in the future.. The easiest and fastest way I found was to switch to PostGIS which has been orders of magnitudes faster and way easier to handle. It didn't end up being all that hard. If you're on a Windows machine, just install PostgreSQL with the PostGIS extension. You'll operate it via PGAdmin. To import geo files to your PostGIS databases you will need the PostGIS importer/exporter which is basically just a window that helps you connect to your Postgres database and upload file to it. There's good youtube videos walking you through the install.Depending what you do it might involves going back in forth between QGIS and PGAdmin. Adding PostGIS layer in QGIS is easy, click on the Elephant icon below the vector and raster "add layer" icon.
Then I basically ran two queries back to back, (running in seconds instead of hours!). The first one buffers each points and dissolve to create clusters.

Create table buffered_cluster_20m (gid serial, geom geometry, area numeric, numID text);
Insert into buffered_cluster_20m(geom, area, numID)
Select (st_dump(st_union(st_buffer(a.geom, 7)))).geom as geometry,
st_area((st_dump(st_union(st_buffer(a.geom, 7)))).geom) as area,
b.numID as numID
From point_samples as a, "donut_20m_cleaned" as b
Where ST_Intersects(a.geom, b.geom)
Group by numID

Problem is, area of cluster is not enough for intended purpose to measure the importance of a cluster (i.e.: 5 points could cover more area than 100 if they are farther apart than each others). So I run a second query that counts the points within the freshly generated areas to get a density instead of just a count:

Create Density_Cluster_20m (bufferID numeric, density numeric);
Insert into Density_Cluster_20m(density, bufferID)
SELECT
(count(a.geom))/(b.area) as density,
count(a.geom) AS counts,
b.gid as bufferID,
b.numid as numID,
b.geom as geom
FROM
"sample_points" AS a,
"buffered_cluster_20m" AS b
WHERE ST_Intersects(a.geom, b.geom)
group by bufferID, numID, b.area, b.geom
order by density desc

And that does it. I end up with a new layer that has areas with a density column that I can import into QGIS and browse through area of interests.


texts: Run a large model by calling it in a python script?

texts:
I've created a rather large model containing a sub-model that calculates the volume of water basins(playas) in a relatively flat area of Texas. The volume is taken from one foot slices of each each basin up to the point of overflow. I would like to use it to batch-process several sets of the input parameters sequentially.

There are 3 input parameters that are as follows: the map extent(XMin YMin XMax YMax) which the user sets around the basin to be processed; the name of a geodatabase to create(String variable); and the name of the basin itself(String) which is referenced within the geoprocessing tasks from an existing feature class. The model works great but I can only run it once.

I would like to be able to use the input as follows:

Move the map view to the extent for the basin and and fill the grid for extent by using the drop down and selecting the "same as display option".
Enter the new name of the geodatabase to be created.
Enter the name of the Playa(basin).
Repeat for all playas that need to be run(Batch).

The problem I run into is that the model will not run sequentially when using it in the batch window. I'm aware that this is a known limitation of the model builder environment and have had many suggestions that what I'm trying to do is better handled by a Python script. So....

Is it possible to write a script that calls the model in which I could batch process from the script using an iterator or some other method to repeat it? If so, how?
I would be okay with entering my parameters into an Excel sheet that could be read by the script. Is this possible? If so, how? Or, is this even necessary?
Would I need to rewrite my model in script form for this to even be possible?

Me:
I have a small knowledge of python scripting, but am well versed in model builder.

I have not done much geoprocessing in python. Mainly small scripts that calculate values and other pretty light tasks.

I've included a picture of the model so you can get an idea of what I'm dealing with.  I know that it won't really be readable. I can provide more in depth examples if needed.


texts: A problem about Correlated-K algorithm of MODTRAN 5

texts:
When I use MODTRAN 5.2.2 to run the test case 'CorrelatedK.tp5', I get the following warning information. How I can solve this issue?

WARNING:  The Correlated-k option was turned off.
           The Correlated-k distributions file DATA/CORK01.BIN

           requires minimally the following assignments in file PARAMS.h:
             MXKSUB =  0        MXGAML = 33
             MXGAMD = 36        MXNUML = 23


texts:
I have downloaded OSGeo4W 64-bit from https://qgis.org/en/site/forusers/download.html and would like to do a 'bleeding edge development build" on my Windows 10. So I select clicked Run and selected Advanced Install, install from internet, set root directory, direct connection, chose the download site, then it gets to the Select Packages stage.

I click Next and immediately get an Installation Status that says "OSGeo4W installation completed successfully" at the top and "Nothing needed to be installed" at the bottom with the Finish button. I have looked around and QGIS isn't installed. What am I doing wrong?


texts: Which version of the Oracle Client do i need to get Oracle support in QGIS 2.6.1?

texts:
I have just installed QGIS 2.6.1 and wish to connect to our Oracle Spatial 11g Database. Am I correct in thinking I will need the Oracle Client installed on my machine for this to work? If so which version?

I'm not sure if it makes any difference but the DB in question is in SDO rather than ESRI's native Binary format.
The Oracle DB version is 11.2.0.3


texts:
I have an error in QGIS that I cannot seem to figure out. I have a zip-code level shapefile from the U.S. Census. I'm merging that with a data file by zip of counts. I want to make a dot density shapefile using the Random Points ftool. However I get a pesky error message as follows:

  "Too many features couldn't be calculated due to conversion error. Please check out message log for more infor"

I've seen the underlying Python code where this error message is created, but I can't seem to find the source.  The message log shows the following

2016-03-21T23:50:23 0   QGIS Ready!
2016-03-22T00:05:17 1   WARNING - fTools: Random Points
2016-03-22T00:05:17 1   The following feature ids should be checked.
2016-03-22T00:05:17 1   Feature id: 33143
2016-03-22T00:05:17 1   Feature id: 33143
....
....
2016-03-22T00:20:04 1   End of features to be checked.

There are many lines in between the dot continuations. So First off, I'm not sure how to check the Feature referenced, as I'm not sure what that feature is.  Without knowing something about the source of the error, I'm not sure how to proceed.

A couple of thoughts. The number of points crammed into a zipcode can be as high as 493. I doubt having too many points in a polygon can create that, but just throwing it out there. Second, most zipcodes are null values. Third, the original Census zipcode polygons were clipped to exclude water along the ocean coasts and the great lakes. Other than that I'm at a loss.


texts: Can GPX Files Contain More Than 1 Track?

texts: Question: can GPX files contain more than 1 track?  I am new to using this format.  Since you can name a track, it suggested to me that I could put multiple tracks in a file.  However, when I tried testing the file out (for exampe, on http://www.mygpsfiles.com/app/), it only displays one of the tracks.


texts: How to use custom offline map in iOS app

texts:
I'm developing an iPhone offline map application. I've illustrator map georeferenced (using mapublisher), have geodatabase (http://goo.gl/isw6IP) and parking data (shp/csv).

Requirements:

The map should Pan, Pinch & zoom.
Search addresses and parking (addresses have alternate names too) on map.
Routing/Navigate to parking/address.
GPS Positioning.

My Illustrator map is a small area of a city and app will work offline without internet connection once installed.

Questions:

In which format I should export the illustrator map, geodatabase &
parking data to be used in a iOS app.
Which iOS offline maps SDK is recommended for above requirements.

Few iOS SDKs like Route-ME, MapBox, Alpstein, Altus, Skobbler, cloudmade, ArcGis iOS SDK, nutiteq etc are using OSM data or their own hosted MAP data but I want to use my own illustrator map and all data will be available local for offline use.

May be above SDKs support what I'm looking for but I don't where to start.


texts:
I want to see comma instead off points in my attribute table.

What happens is: When I import TXT or CSV there´s an option to use comma as decimal delimiter but after the importation I only see points as decimal delimiter.

What do I need to do?

Source CSV

Before

After


texts: Exposing Active Directory Groups on ArcGIS Server 10.2.2

texts: Is there a way to have ArcGIS Server 10.2.2 expose Active Directory Groups and users instead of just users?  I have IIS and the Server set up for Windows Authentication (IIS) and Web Authentication (ArcGIS server).  I'm looking to secure services whereby only members of certain Active Directory groups we have set up in the organization can access.  This would be preferred to the only option I have seen so far - creating a role wmanually selecting each user


texts:
Our road is sometimes (cannot reproduce now) marked with red stripes on Google maps. Not sure what that means, but at least it does not navigate to this road.

This means that Google Maps leads people to the nearest freeway, and indicates they have arrived when they are on the freeway adjacent to our house. This is both annoying and dangerous in cases where people try to stop or such.

I have tried submitting map changes many times, during the previous 5 years or so, but nothing is happening. I am unable to find somebody to contact in this regard, and very frustrated that this can apparently not be fixed.

Any ideas what I could do to fix this?


texts: Splitting A polygon into multiple polygon by multiple line strings in Leaflet and turf.js

texts:
I have need to split a polygon into multiple polygons based on the line strings drawn on the polygon. I have come across another post that splits a polygon into two based one line string. If a new line is drawn then the existing lines are ignored and the polygon is split into two based on the new line. What I want to achieve is 'split a polygon into multiple polygons based on multiple line.

 var osmUrl = 'http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png';
  var osmAttrib = '&copy; OpenStreetMap contributors';
  var osm = L.tileLayer(osmUrl, { maxZoom: 18, attribution: osmAttrib });
  var drawnItems = L.featureGroup();

  var map = new L.Map('map', { center: new L.LatLng(51.505, -0.04), zoom: 13 });

  osm.addTo(map);
  drawnItems.addTo(map);

  map.addControl(new L.Control.Draw({
    draw: {
      marker: false,
      circle: false,
      circlemarker: false,
      rectangle: false,
      polygon: {
        allowIntersection: false,
        showArea: true
      }
    }
  }));

  function cutPolygon(polygon, line, direction, id) {
    var i = -1;
    var j;
    var polyCoords = [];
    var retVal = null;

    if ((polygon.type != 'Polygon') || (line.type != 'LineString')) return retVal;
    if (line.coordinates.length != 2) return retVal;

    var intersectPoints = turf.lineIntersect(polygon, line);
    var nPoints = intersectPoints.features.length;
    if ((nPoints == 0) || ((nPoints % 2) != 0)) return retVal;

    var offsetLine = turf.lineOffset(line, (0.01 * direction), {units: 'kilometers'});
    var thickLineCorners = turf.featureCollection([line, offsetLine]);
    var thickLinePolygon = turf.convex(turf.explode(thickLineCorners));

    var clipped = turf.difference(polygon, thickLinePolygon);

    for (j = 0; j < clipped.geometry.coordinates.length; j++) {
      var polyg = turf.polygon(clipped.geometry.coordinates[j]);
      var overlap = turf.lineOverlap(polyg, line, {tolerance: 0.005});
      if (overlap.features.length > 0) {
        polyCoords[++i] = turf.coordAll(polyg);
      };
    };

    if (i == 0)
      retVal = turf.polygon(polyCoords, {id: id});
    else if (i > 0) {
      retVal = turf.multiPolygon([polyCoords], {id: id});
    }

    return retVal;
  };

  var polygons = [];
  var layers = [];

  map.on(L.Draw.Event.CREATED, function (event) {
    var layer = event.layer;
    drawnItems.addLayer(layer);

    var geojson = layer.toGeoJSON();
    var geom = turf.getGeom(geojson);

    if (geom.type == 'Polygon')
      polygons.push(geom);
    else if (geom.type == 'LineString') {
      var line = geom;
      layers.forEach(function (layer, index) {
        layer.remove();
      });
      layers = [];
      polygons.forEach(function (polygon, index) {
        var layer;
        var upperCut = cutPolygon(polygon, line, 1, 'upper');
        if (upperCut != null) {
          layer = L.geoJSON(upperCut, {
            style: function(feature) {
              return {color: 'green' };
            }
          }).addTo(map);
          layers.push(layer);
        };
        var lowerCut = cutPolygon(polygon, line, -1, 'lower');
        if (lowerCut != null) {
          layer = L.geoJSON(lowerCut, {
            style: function(feature) {
              return {color: 'yellow' };
            }
          }).addTo(map);
          layers.push(layer);
        };
      });
    };
  });

I want the polygon split into 12 polygons. now it's just two(green & yellow).

I got this from this here


texts: Reasons of displacement between same photo taken from two different satellites

texts: So I am interested what could be the reasons for displacement when comparing the same photo taken from two different sattelites (ASTER, LANDSAT). I need this for my study so I would be thankful for getting an answer.


texts: QGIS - Join polygon layer attributes to line segments

texts:
I've got a hexagon layer with data that I want to have joined to the available road network for the area so that the road segments have the attributes from the hexagon layer.
Each segment of the road layer is comprised by the intersection of two roads, there are a total of 304 features. The features that go through multiple polygons should account for the average value of these polygons. Say that H_1 = 0.4 and H_4 = 0.2 (whereas H -> hex polygon bin) thus, (H_1+H_4)/n (where n the number of polygons that intersect the line layer).
How can I do that?


texts:
I am using QGIS 3.10.0

I am attempting to join a spatial vector layer to a non-spatial layer (CSV).

Each polygon on the spatial layer has a unique id attribute.  The CSV also has this id attribute but it is not unique for each record.  So for each polygon on the spatial layer the are 1-to-N associated records in the CSV.

I was hoping that when you select an object on the spatial layer it will show all the associated records from the CSV.

I initially used the layer join option but that only shows one record from the CSV.

I then read about virtual layers.  I created a virtual layer of both existing layers and this does what I want... however this is so slow it is unusable!  After creating the layer it takes 10 minutes for QGIS to become responsive.  If I do anything with that layer (select an object or just zoom!) I have to wait another 5 to 10 minutes!  This is on 64-bit machine core i7 16gb ram.

The spatial layer has 7988 records and the csv has 8341 records.

Does anyone have any other options or know how to speed the virtual layer rendering up?


texts: PostGIS: Simple query selecting nodes of polygon

texts:
Asked also at stackoverflow

INPUT:

Regular OSM buildings .shp table
PostGIS Point location
RADIUS (around above location), lets say 500m.

RESULT:

I want to get back all corners from all buildings within this radius from location, that ("WHERE") they distance from location is less than DISTANCE (let's say 200m).

I'm struggling to get result as points accepting WHERE clause from polygons table.

EDIT:

I started something like this:

SELECT ST_AsText((ST_DumpPoints(geom)).geom)
FROM buildings
WHERE ST_DWithin((ST_DumpPoints(geom).geom), ST_GeomFromText('POINT(5 45)',4326),  0.002);


texts: Extracting list of coordinates of all pixels in neighborhood of given pixel in R

texts: I have a Raster layer. I want to get coordinates of all pixels in neighborhood (for example 7*7) of given pixel. Is there any function available to do this in R?


texts: PyQGIS changeAttributeValue not saving modified records

texts:
I've been trying to save some changes to the a dataset using PyQGIS but I can't see any changes even after closing/reloading the layer attribute table.
I have tried a few different options but none of them work.
import re
layer = iface.activeLayer()
selection = layer.getFeatures() #all features
#selection = layer.selectedFeatures() #only selected features
layer_provider=layer.dataProvider()
layer.startEditing()
for f in selection:
   text=f['RoadMntnc'] #get the unique road id to process.
   print (text)
   stext=re.split('_+', text) #split the road id into it's components.
   print(stext)
   rtype=stext[-1:] # the type of road
   newtext=' '.join(stext[:-1]) # the road name without the road type
   print (newtext)
   attr_value={2:newtext} # the column to replace and the value to use
   #layer_provider.changeAttributeValues({1:attr_value}) #attempts to write the value
   #layer.changeAttributeValue(f.id(), RoadName, 30)
   layer.changeAttributeValue(f.id(), 2, newtext)
layer.commitChanges()
layer.updateFields()

Dataset looks like

I want to eventually autofill the RoadName, RoadType and other fields from just the unique RoadMntnc name so that there's less work/less data entry issues when adding roads to the dataset.


texts:
Is it possible to fill no data around bridges without pixels blurriness (see picture). I use Elevation Void fill function, also tried in qgis No data filling and Grassgis r.fillnulls. I know, that I can just assing all Nodata values to specific value, but I need to get smooth raster.

before

after


texts: GeoNetwork + Mapproxy customized endpoints?

texts:
I'd like to build a service that allow each user of a GeoNetwork portal or sub portal (https://geonetwork-opensource.org/manuals/trunk/eng/users/administrator-guide/configuring-the-catalog/portal-configuration.html) to share maps.

By default, we can share maps with OWS context, But the fact is that client applications must be modified to handle the OWS context.

Do you know a plugin or approach that can enable the publication of a per-portal or per-subportal decated sets of endpoints (WMS/WMTS/WFS)?

Perhaps by associating MapCache/Mapproxy with GeoNetwork?

Therefore, this will allow immediate compatibility with client applications, allowing them to use simple GetCapabilities requests on standard endpoints.


texts: Generic Writer Fanout of file geodatabase

texts:
I would like to fanout file geodatabase on an specific attribute by using generic writer with Dynamic schema. I know how to do this when not using Generic Writer, so I assumed that principle will be the same. What I get now are fanout folders comprising of fanout fgdbs, rather than only fanout fgdbs. I.e. in the Fanout Directory are folders (FanoutPrefix_AttributeTOFanoutOn) and in those folders are fgdbs (BaseFilename_AttributeToFanoutOn).

So, I do not need to have fanout fgdb to be placed in the additional folders created when doing fanout, only to be placed in the Fanout Directory as defined in Navigator.


texts:
Azure SQL has a hemisphere limitation where it errors on any geography having an envelope angle greater than 90 degrees.

My coworker wants to set up some detection in the web server to catch such instances before they hit the database. How can envelope angle be calculated from a series of coordinates that define a polygon or multipolygon?


texts: Why provider.nextFeature(feat) does not iterate over all layer features?

texts:
I have a point layer with over 1700 features, but iterating this way:

    i = 0

    while provider.nextFeature(feat):
        i =i + 1
        print '{0} de {1}: {2}'.format(i, provider.featureCount(), feat.id())

...I get only a few of them and no error arises.

Could it be possible to make a list of feature ids and then iterating over each QgsFeature by their id()?


texts: GeoServer: is grades to degrees transformation possible through SLD file?

texts:
I work on GeoServer 2.6.2 and PosgreSQL 9.2 / PostGIS 2.0. I have declared a layer in GeoServer which is linked to a PostgreSQL database table. In this table, I have a column (obj_orientation) storing orientation values in grades (360° = 400 gr.). Is it possible to transform these values to degrees inside to SLD stylesheet?

I have tried these but the layer doesn't display at all:

<Rotation>
   <ogc:PropertyName>obj_orientation*0.9</ogc:PropertyName>
</Rotation>

<Rotation>${obj_orientation*0.9}</Rotation>

<Rotation>${obj_orientation}*0.9</Rotation>

Is there a way to do so or do I have to update my column in the database?


texts: Seeking open source web client alternative to Geocortex Essentials?

texts:
What is best the opensource webclient alternative To Geocortex essentials?
For someone  without much coding experience.


texts: Selecting random points from specific attribute in yield map using QGIS

texts:
I am using QGIS.

A farmer has sent me the shapefiles from his John Deere combine. They were three shapefiles that I have merged into one layer. The shapefile is from one field with a treated tramline and a control tramline. I have georeferenced an image that shows where the treated points are and where the control points are.

Following that, I have selected the features from both tramlines individually and created a separate layer for each one.

I would like to select random points from each layer and do some spatial statistics.

I have tried using the minimum bounding geometry to create a convex hull around the layer but the attributes don't get copied when I do so which means that I am not able to use the random points in layer function.

The attribute table consists of 10 variables - the one I am interested in is wet mass.

The wet mass attribute is not an integer number.

How do I select random points from a points layer that in its attribute table has a number of variables but I only want to use one that the type of value is REAL

This is what I see when I open the attributes table of the created boundary layer


texts: I'm trying to understand what is the point of GeometryField.  The GeoDjango docs don't say much and neither do questions on Stack.  What is special about it?  When would you use it over PolygonField/PointField/etc.?


texts: How to store large time series of raster images

texts:
I've got a large time-series of rasters, derived from satellite images. Each raster is about 2000 x 2000 pixels, and I've got about 10,000 individual images. I'm trying to work out what the best way to store these is.

Of course, the best way to store them will depend on what I want to do with them - and in this case I want to be able to look at a time series for individual pixels, time-series averaged over specific areas (eg. vector administrative area data overlain on the raster), produce average images over various timescales, and view individual images.

The approach I've used before for shorter time series has just been to stack all of the images into a multi-band GeoTIFF - but I'm not even sure if a GeoTIFF will allow you to have 10,000 bands, or how well it will perform. An alternative would be to store them in thousands of separate files, but that would probably require a lot of custom programming work to be able to do the analysis.

I'm sure I can't be the first person to have this problem - what approaches are recommended for dealing with this sort of volume of time-series data?

Ideally I'd prefer solutions that work nicely with Python and don't require ArcGIS - but I'd be interested in any sensible ways forward.


texts: ArcMap event to detect when a toolbar open/close with ArcObjects?

texts:
I have this code to know if a toolbar is displaying:

IApplication m_aplication = ArcMap.Application;

ICommandBar versionToolbar = GetToolbarByName(m_aplication, "esriArcMapUI.VersionToolbar");

if (versionToolbar.IsVisible()){
    versionToolbar.Dock(esriDockFlags.esriDockHide);
}

I can hide a toolbar, but I need to know when a toolbar is added to the desktop view, for example when the user check a toolbar in customize option:

I have this extension class:

public sealed class eManager : ESRI.ArcGIS.Desktop.AddIns.Extension
    {
       protected override void OnStartup()
       {
            MessageBox.Show("on Start up eManager");
       }
    }

And this config extension:

but the event doesn't work


texts: Is there a Python Lib for requesting WMS/WFS and saving as image/PDF?

texts:
i am wondering whether there's Python open-source GIS lib which has APIs to support call WMS/WFS from another GIS server (e.g., GeoServer) and then save the response data(WMS Basemap and WFS layer) as pictures.

any recommendations?

thanks for  any inputs!

UPDATE:

what I am trying to do is a Map Printing service, by using OpenLayers as the front-end and Django as the server; Client user set the extent and layers and then send the print request (which refers to the parameters, i.e., map extent, names of layers) to server, then server takes over this request and call WMS/WFS again by using request parameters, save the response as PDF, export this PDF link to client.

The difficult part is that how the server call WMS/WFS and combine/overlay these responses together (i.e., put these map/layers together, since WMS is usually the base map, WFS points to the feature layers), finally save this combined object as Image.

in current answers, urllib seems a good one, but i am not sure how to combine these responses (WMS, WFSs) together; OWSLib also seems another good option, but it indicates it's a client programming tool, I am a little confused that whether it's appropriate for my use...

any other further inputs???

appreciate!


texts: Why is the SQL included in the URL of CartoDB visualizations?

texts:
I noticed that the URL of CartoDB visualizations include the parameter sql= with the SQL statement for that visualization:

http://lifewatch-inbo.cartodb.com/viz/b43fa5a6-7dd0-11e3-ba2a-5f9e077d1162/embed_map?title=true&description=true&search=false&shareable=true&cartodb_logo=true&layer_selector=false&legends=true&scrollwheel=true&sublayer_options=1&sql=SELECT%20%0AST_MakeLine(the_geom_webmercator%20ORDER%20BY%20date_time%20ASC)%20AS%20the_geom_webmercator,%20%0ACASE%0AWHEN%20date_time%20%3E%3D%20'2013-05-28'%20AND%20date_time%20%3C%20'2013-06-25'%20THEN%20'1.%20Egg%20incubation'%0AWHEN%20date_time%20%3E%3D%20'2013-06-25'%20AND%20date_time%20%3C%20'2013-07-24'%20THEN%20'2.%20Chick%20rearing'%0AELSE%20'3.%20Postfledging'%0AEND%20as%20behaviour%0AFROM%20three_gulls%0AWHERE%20bird_name%20%3D%20'Eric'%0AAND%20outlier%20IS%20NULL%0AGROUP%20BY%20behaviour%0AORDER%20BY%20behaviour&sw_lat=50.790527850939384&sw_lon=2.2449874877929688&ne_lat=51.76975203254081&ne_lon=5.587577819824219

However, these parameter seems optional, as the visualization works fine without the parameter:

http://lifewatch-inbo.cartodb.com/viz/b43fa5a6-7dd0-11e3-ba2a-5f9e077d1162/embed_map?title=true&description=true&search=false&shareable=true&cartodb_logo=true&layer_selector=false&legends=true&scrollwheel=true&sublayer_options=1&sw_lat=50.790527850939384&sw_lon=2.2449874877929688&ne_lat=51.76975203254081&ne_lon=5.587577819824219

Changing the SQL on CartoDB will update the existing visualization and render the value of SQL= outdated, so I was wondering why it is included at all?


texts: Reprojecting a raster results in larger raster values than the original raster values

texts: When I reproject a raster in ArcMap using the project raster function (details in the image above) to change its X,Y coordinate system from North_America_1983_Transvers_Mercator to NAD_1983_Idaho_TM, the values in the result raster tend to increase than the original raster values. How can I fix this?


texts: Viewing Revit file using ArcGIS Javascript API

texts: I have a Revit file with GPS information in it. Is there any provision to view that Revit model in ArcGIS Web directly without using ArcGIS desktop    by publishing it as a web view and then opening that web view.


texts: Is there an open source directions webclient for pgRouting like the directions API from Google?

texts:
I am looking for an open source webclient with (Google) directions functionality (like draggable lines to create extra waypoints) based on an pgrouting database (with an osm2pgrouting topology).

I noticed that there are some examples and demos on the pgrouting site (some are not working anymore), but they are all specific things people have done with pgrouting. They seem to have all created their own specific front end. I am looking for a more generic directions/routing client based on a pgrouting topology.


texts: Selecting features within polygon from another layer in python?

texts:
Someone asked a nearly identical question about selecting features in a layer by location here Selecting features within polygon from another layer using R or Python?, but the answers are not for general gdal/ogr python users. I'm wondering if anyone has a pythonic solution to this problem that does not require qgis modules. In particular, I'd like to be able to specify the options qgis offers (include input features that touch and/or overlap/cross and/or are completely within and/or are selected).

I've been trying all day to remove polygons from a shapefile that share a border with (but are contained completely within) a polygon in another shapefile. This task takes <5min in QGIS, but after 5+ hours I still haven't found a way to "select by location" and specify "options" using just ogr (I've tried using layer.SetSpatialFilter(polygon_geometry), but the docs specify that this only works with GetNextFeature(), which I am not using to loop through my features, so I was not terribly surprised when it did not produce the desired results. Also, I suspect that features sharing only a border with the selection layer are included rather than excluded). I've also tried ogr2ogr

ogr2ogr -f 'ESRI Shapefile' -dialect sqlite -sql "SELECT * FROM bergs WHERE ST_Intersects(bergs.geometry, AOI_buffered100.geometry)" -overwrite bergs_testborderrem.shp bergs.shp

but cannot get a working SQL statement (I keep getting no such column errors, but all the ogr-sql examples I've found use .geometry, so I have no idea how to query on the geometry otherwise).

Does anyone have any suggestions?