texts: GEE JavaScript Expecting FeatureCollection but receiving Feature

texts:
I am working on calculating NDVI values for both SEN2 and SEN3 from the same time frame and storing the results on my Google Drive.
I have a feature collection called "region" which is an assemble of 4 fields (polygons) and is imported from my assets.

I run this code to compute and save NDVI values for each field and for both SEN2 and SEN3 but I am unable to save the outputs as the task cannot be complete due to this error:

Error: Collection.geometry, argument 'collection': Invalid type. Expected type: FeatureCollection. Actual type: Feature.

This is the code:

//define all AOI
var field1 = ee.FeatureCollection(region.filterMetadata("system:index","equals","00000000000000000000").first());
var field2 = ee.FeatureCollection(region.filterMetadata("system:index","equals","00000000000000000001").first());
var field3 = ee.FeatureCollection(region.filterMetadata("system:index","equals","00000000000000000002").first());
var field4 = ee.FeatureCollection(region.filterMetadata("system:index","equals","00000000000000000003").first());
var fields = [field1, field2, field3, field4];

//Get the ndvi for the specific field and time
function get_ndvi_sen2(region, startdate, enddate){
  region = ee.FeatureCollection(region)
  var ndvi_s2 = ee.ImageCollection('COPERNICUS/S2_SR')
    .filterDate(startdate, enddate)
    .map(function (image) {
      return image
        .clipToCollection(region)
        .normalizedDifference(['B8', 'B4'])
        .updateMask(image.select('QA60').not());
    })
    .median()
    .rename('ndvi_s2');

  var ndvi_s2_Features = ee.FeatureCollection(ee.List(ndvi_s2
  .reduceRegion({
    reducer: ee.Reducer.toList(),
    geometry: region,
    scale: 10
  })
  .get('ndvi_s2'))
  .map(function (ndvi_s2) {
    return ee.Feature(null, {ndvi_s2: ndvi_s2});
  })
  );
  return ndvi_s2_Features;
}
function get_ndvi_sen3(region, startdate, enddate){
  var SEN3 = ee.ImageCollection('COPERNICUS/S3/OLCI')
                    .filterDate(startdate, enddate);

  // Select bands for NDVI and apply band-specific scale factors.
  var SEN_corr = SEN3.select(['Oa17_radiance', 'Oa08_radiance'])
                .median()
                // Convert to radiance units.
                .divide(ee.Image([0.00493004, 0.00876539]));
  region = ee.FeatureCollection(region)
  var ndviSEN3 = SEN_corr
        .clipToCollection(region)
        .normalizedDifference(['Oa17_radiance', 'Oa08_radiance'])
        .rename('ndviSEN3');

  var ndviFeatures = ee.FeatureCollection(ee.List(ndviSEN3
    .reduceRegion({
      reducer: ee.Reducer.toList(),
      geometry: region,
      scale: 10
    })
    .get('ndviSEN3'))
    .map(function (ndviSEN3) {
      return ee.Feature(null, {ndviSEN3: ndviSEN3});
    })
  );
  return ndviFeatures;
}
function write_to_drive(col, disc, selc, folder){
  Export.table.toDrive({
    collection: col,
    description: disc,
    selectors: selc,
    folder: folder
  })
}

var start = '2019-03-1'
var end = '2019-04-01'
var sen2_disc, sen2_selc = 'ndvi_s2';
var sen3_disc, sen3_selc = 'ndviSEN3';
var sen2_folder = 'GOOD/SEN2'
var sen3_folder = 'GOOD/SEN3'

for (var i in fields){
  var sen2_loop = get_ndvi_sen2(fields[i], start, end);
  var sen3_loop = get_ndvi_sen3(fields[i], start, end);
  write_to_drive(sen2_loop,sen2_disc,sen2_disc,sen2_folder)
  write_to_drive(sen3_loop,sen3_disc,sen3_disc,sen3_folder)
}


texts: Where can I host GeoServer that will allow multiple users to edit on ArcMap?

texts: We have several spatial data editors spread across the US, each having to access the same data. What we need is a multi-user spatial database that can be served up through GeoServer, hosted on the cloud, and accessed/edited within ArcMap. Does this exist?


texts: Accessing array of values(height, intensity) associated with raster from ArcObjects without looping?

texts:
I am trying to access the array of point values(mainly height and intensity) from a raster through C# without looping through values. I am looking for a direct pointer to the array of values inside of the Raster. I have already tried using methods such as IRaster2.getVal and ConversionOp.RasterDataToPointFeature data but they all take way to longand we should just be able to access the Array of points directly. When using PixelBlock3.get_PixelData we recieve an array filled with zeros. This could be because I am not casting the Array correctly.

Here is the source for the function I am trying to write:

    public double[,] RasterizeLASDByEnvelope(string Name, IGeometry geometry, IRaster pRaster)
    {
        Debug.WriteLine("[RasterizeLASD]");
        Stopwatch sw = new Stopwatch();
        Debug.WriteLine("RasterizeLASD: Start Process... ");
        sw.Start();
        Debug.WriteLine("RasterizeLASD: Defining Geoprocessor.");
        IGeoProcessor2 gp = new GeoProcessorClass();
        gp.OverwriteOutput = true;
        IGeoProcessorResult2 result = new GeoProcessorResultClass();
        IVariantArray parameters = new VarArrayClass();

        object sev = null;
        try
        {

            IRasterProps pRasterProps = pRaster as IRasterProps;

            Debug.WriteLine("RasterizeLASD: Initializing Point Value
            Lists...");
            int width = pRasterProps.Width;
            int height = pRasterProps.Height;
            Stopwatch cpsw = new Stopwatch();
            Debug.WriteLine("RasterizeLASD: Start Point Transfer...");
            cpsw.Start();
            IConversionOp conversionOp = new
            RasterConversionOpClass();
            IFieldChecker fieldChecker = new FieldCheckerClass();
            string name = "pnt" + Program.gdbAdapter.DateTimeStamp();
            string finalname;
            fieldChecker.ValidateTableName(name, out finalname);
            IGeoDataset featClass = conversionOp.RasterDataToPointFeatureData((IGeoDataset)pRaster, Program.gdbAdapter.getWorkspace(), finalname);
            IFeatureCursor fc = ((IFeatureClass)featClass).Search(null, false);
            IFeature feature = fc.NextFeature();
            double[,] solution = new double[width, height];

            for (int y = 0; y < height; y++)
            {
                for (int x = 0; x < width; x++)
                {
                    object pntval = feature.Value[3];
                    double value = (double)pntval;
                    solution[x, y] = value;
                    feature = fc.NextFeature();
                }
            }
            Debug.WriteLine("RasterizeLASD: End Point Tranfer... Transfer Time: " + cpsw.Elapsed);
            sw.Stop();
            Debug.WriteLine("RasterizeLASD: End Process... Process Time: " + sw.Elapsed);
            return solution;
        }
        catch (Exception ex)
        {
            // Print a generic exception message
            Debug.WriteLine(ex.Message);
        }
        return null;
    }


texts:
I have a sfc-object consisting of just points. Unfortunately, when I read in the shapefile, it created points of type XYZ. Now I removed some points and can't figure out how to export the remaining points as shapefile. I get the error:

GDAL Error 6: Geometry type of `3D Point' not supported in shapefiles

As I don't really need the third dimension So I thought I'd just remove it. But I don't know how. This is how the sfc-object looks like at the moment.

First 5 geometries:
POINT Z (397864 6302484 0)
POINT Z (386440.2 6299046 0)
POINT Z (390490.3 6296620 0)
POINT Z (391614.7 6293960 0)
POINT Z (377803.1 6276244 0)


texts: Create a map survey and download the result

texts: I have a shapefile of our state boundary. I want people to be able to drop a pin and add their names to the pin as an attribute. From there I would like to download all the pins as a shapefile and import it into gis software. Is anyone aware of a website or app that could accomplish this. I do not need a basemap or imagery as I want the users to drop the pins on where they think the location is at without any context. I'm not familiar with Leaflet, but could that possibly be a solution?


texts: Adding weather information to QGIS atlas composer

texts: I'd like to add a tabular weather forecast with climatic variables like temperature, wind, and chance of precipitation to my QGIS atlas composer.  Does anyone have suggestions on how I can do this?


texts: How to design markers - raster or vector?

texts:
I have several layers that are point geometries. I am defining them custom markers symbols and such symbols have colors (red, yellow and green) to represent values.

Since there are several of them in an area, I started with marker clusters, but quickly realized that clusters are good only if all markers mean the same (like redfin home markers).

I would like to stick to markers but to avoid visual clutter - I would like to resize them based on the zoom level. If I design these markers in SVG, I could achieve this - as its easy to size SVG than raster formats, similar to this link http://mapbox.com/tilemill/docs/guides/advanced-map-design/#images_as_icons

Please send me any examples or sample sites you have developed or come across so that I can see the visual appeal of having several custom markers in maps.

Thanks,
Ramesh


texts: Google maps satellite view Pricing and limitations

texts:
This script at openstreetmap shows how one can use google's Satellite view as a layer on top of others.

I noticed that to do so, one does not need to supply any API key. Does that mean that there are no limitation to the number of satellite views requests one can do?

In my case, I would like to use leaflet to display layers from both mapbox  and google's Satellite view. But the cost of using the later is unclear!

I hope I made my point clear enough!

Thank you.


texts: Swapping raster bands on disk efficiently using R

texts:
I have two raster files, each with RGB bands and an alpha layer. I would simply like to swap one of the raster's alpha bands for another, but to do so as efficiently as possible (these are very big files and there are a lot of them)

I started by trying this:

library(raster)

ras1 <- brick("raster1.tif")
ras2 <- brick("raster2.tif")

ras1[[4]] <- ras2[[4]]

writeRaster(ras1, "newRaster.tif")

This, I believe, ends up rewriting the entire four layers, even if only one of the bands has changed.

I then thought I had found the perfect function with raster::update, but this turns out to be extremely slow

update(ras1, getValues(ras2[[4]]), 1, 4)

I am assuming it is because using getValues reads a bunch of data into memory. Can anyone suggest a more efficient way of doing this?


texts: Generate MVT with POSTGIS and feature tolerance

texts:
I have a database which has three geometry tables 'tbl_point, tbl_line, tbl_polygon' with millions of features.
Each table has two geometry columns named geom & m_geom for geometries and multi-geometries respectively.
I am generating MVT using POSTGIS with following query.
WITH
   bounds AS (SELECT ST_Segmentize(ST_MakeEnvelope(8247861.100083642, 3546678.1124321707,
      8252753.069893893, 3551570.082242422, 3857),1222.9924525627866) AS geom,
      ST_Segmentize(ST_MakeEnvelope(8247861.100083642, 3546678.1124321707, 8252753.069893893,
                    3551570.082242422, 3857),1222.9924525627866)::box2d AS b2d
     ),
  f_300 AS (SELECT ST_Simplify(ST_AsMVTGeom(ST_Transform(t.geom, 3857), bounds.b2d, 4096, 128, true), 20)
            AS geom, t.id FROM tbl_line t, bounds
            WHERE ST_Intersects(t.geom, ST_Transform(bounds.geom, 4326))
            AND t.shp_id = 300 AND t.deleted = FALSE
            UNION
            SELECT ST_Simplify(ST_AsMVTGeom(ST_Transform(t.m_geom, 3857), bounds.b2d, 4096, 128, true), 20)
            AS geom, t.id FROM tbl_line t, bounds
            WHERE ST_Intersects(t.m_geom, ST_Transform(bounds.geom, 4326))
            AND t.shp_id = 300 AND t.deleted = FALSE
           ),
 f_299 AS (
           ......
           AS geom, t.id FROM tbl_point t, bounds
           ......
          ),
f_307 AS ( ......
           AS geom, t.id FROM tbl_polygon t, bounds
           ......
         )

SELECT ST_AsMVT(f_300.*, '300') FROM f_300
UNION
SELECT ST_AsMVT(f_299.*, '299') FROM f_299
UNION
SELECT ST_AsMVT(f_307.*, '307') FROM f_307

The query exludes small polygon using ST_Simplify tolerance. But for points and lines its returning all the features which becomes very congested while displaying.
Is there any way to show only some of the features or group features or cluster features.
Also, any recommendation on improving the query is highly appericiated.


texts: pandas.to_csv to write arcpy.getrasterproperties returns

texts:
I'm stumped on how I can use pandas.to_csv for this application. I'm trying to first write a CSV then append the CSV within the for loop so that every iteration appends a new row.

import arcpy, numpy, pandas
from arcpy import env
env.workspace = r"C:\VMshared\small_example_valley5\SDepth2"
headers = ["ROWCOUNT", "COLUMNCOUNT", "MAX", "MIN"]
pandas.to_csv(r"C:\VMshared\small_example_valley5\output_tabs\CalcPropOut.csv", header = headers, mode = "w")

for file in arcpy.ListRasters():
    ROWCOUNT = arcpy.GetRasterProperties_management(file,"ROWCOUNT")
    COLUMNCOUNT = arcpy.GetRasterProperties_management(file,"COLUMNCOUNT")
    MAX = arcpy.GetRasterProperties_management(file,"MAXIMUM")
    MIN = arcpy.GetRasterProperties_management(file,"MINIMUM")
    pandas.to_csv(r"C:\VMshared\small_example_valley5\output_tabs\CalcPropOut.csv", sep=",", na_rep='-9999', columns=headers, mode="a")


texts: How to embed ArcMap graph into winforms dialog using .NET/C# and ArcObjects?

texts:
I am creating an addin for ArcMap in csharp and want to add a graph to a dialog.

The native arcgis graph control doesn't seem to been available from the visual studio toolbox.

I can't figure out how to do this from code.

How can i do this?


texts: How to plot a point on a static Google map [png]?

texts:
http://maps.googleapis.com/maps/api/staticmap?center=Brooklyn+Bridge,New+York,NY&zoom=13&size=600x300&maptype=roadmap&markers=color:blue%7Clabel:S%7C40.702147,-74.015794&markers=color:green%7Clabel:G%7C40.711614,-74.012318&markers=color:red%7Ccolor:red%7Clabel:C%7C40.718217,-73.998284&sensor=false

shows a static map which is a png file.

If I do have a coordinate which lies between the coordinates of this static map file, then how should I plot that point on this map (which is basically a png)?

EDIT:

I would be doing it offline. Can't use the maps API. Will save that map as a PNG, then I have to plot the points.


texts: Calculating land sea area ratio in a portion of raster map in QGIS

texts:
I wanted to know how can I calculate land/sea area ratio in a portion of raster map (attached in the mail)
Till now i had this raster map that i converted in to vector by the option raster>>conversion>>polygonise then i install the plugin CATOGRAM to draw grids. However i get no output.

. .


texts: Create polygons from points representing the corner vertices (ArcGIS)?

texts:
I am on ArcGIS 10.3.1. I have a point feature class that represents the corner vertices of the polygons. They are located in a regular grid. The labels represent the OBJECTID. The points can be digitized in a random order with no sequence whatsoever. There are no attributes defined.

What I need to do is to create individual polygons (essentially connecting points with certain nearly located points). If I will be able to connect them as on the figure below, I would create polygons from those lines (ie polygons (17,18,21,22);(18,19,22,23) and so forth). The sequence of vertices in the polygons doesn't matter.

Since those polygons can vary in the shape (can be squares or rectangles with varying side size), I have a hard time finding an algorithm to connect only certain points (in human terms, only horizontally or vertically). Relying on the lat/long is not an option since the points may not lay exactly on the same coordinate lines (ie, may be shifted as below). However, all points in "rows" and "columns" lay on the imaginary line that would go through the edge points of the grid (ie, points 22,23,29 lay on the line that go from 21 to 24).

I am trying to automate the process fully, but if certain steps that don't require much time would be needed to get done manually, it's acceptable. In other words, I am looking for any way to speed up this process.

I use arcpy and Python techniques for finding the 5 nearest points for each points and connect them by creating a new line segment. As you can see, some of the lines created go as diagonals within the rectangles created and should be cleaned up.

Things can get harder when the grid has varying the distance between columns/rows such as on this picture (point 22 gets connected not only to 28, but also to point 27 - outside the rectangle created).

I have looked at Parcel Editor in ArcGIS - couldn't find any option that would let me to create polygons automatically. I have a feeling that finding out the sequence of filtering the lines & resultatnt polygons and re-iterating the output I could get sufficiently good results.

Does anyone have any good tips on how to approach this problem to create polygons from the points representing the corners?


texts: Bing maps get Distance between two places

texts:
I have referred the sample http://help.arcgis.com/en/webapi/silverlight/samples/start.htm#BingRouting and developed a Bing maps routing tool on my Silverlight GIS website.. now my client wants me to also show the distance between two places whos direction is shown..

I took a look at the Routing class but failed to find anything to get the distance.. I used the ‘haversine’ formula, but it does not give me the exact distance along the route..

or

Can I calculate the length of the geometry/polyline drawn on the map.. Will this result be accurate? Please guide me so that I can get the distance between my start and destination..


texts: Exporting Feature Classes to Shapefiles in individual folders of same name using ArcPy?

texts:
I have a Geodatabase full of featureclasses.
I want to perform the following:

Copy out each individual featureclass as an output individual shapefile (and have the same name).
Create individual Output Folders (of the same output shapefile name).
Copy each output shapefile into its own unique folder (of the same name).

Currently, I can perform task 1 and 2 above, however, I am having trouble with task 3 above.  My results so far is that I can create the output shapefiles and folders (of the same name), however, the shapefile outputs do not sit in the pertaining folders. i.e. My result = a bunch of shapefiles and a bunch of folders.  I need to (i) either have the shapefiles created and automatically placed into the pertaining folder (based on the same name of the shapefile), or (ii) After creating the shapefiles and folders, move or copy the output shapefiles into each pertaining folder (based on the same name as the shapefile).

Therefore, for each of my featureclasses, my end result would be 1 seperate folder housing 1 shapefile (same name for both folder and shapefile).

I am fairly new to Python.

# import system modules
import arcpy
import os
from arcpy import env

# Set environment settings
# Note: The INPUT Geodatabase below houses 5 x featureclasses
env.workspace = r"D:\INPUT.mdb"
arcpy.env.overwrite = True

# Set local vairables
# Note: Ensure that you have this empty folder already created within FileManager etc.
outWorkspace = "D:/OutputFoldersAndSHP"

print "Commencing"

# Use the ListFeatureClasses function to generate a list of featureclasses in the env.workspace further above
fcs = arcpy.ListFeatureClasses()

for fc in fcs:
    # PHASE 1: CREATE A FOLDER BASED ON INDIVIDUAL FC NAMES
    # Note: (It simply uses the fc name as the output folder name!)
    arcpy.CreateFolder_management(outWorkspace, fc)

    # PHASE 2: COPY FEATURECLASSES TO SHAPEFILES
    # Determine the new output shp path and name
    outSHP = os.path.join(outWorkspace, fc)
    arcpy.CopyFeatures_management(fc, outSHP)

print "Completed"


texts: Running Maximum Likelihood Classification?

texts:
I'm trying to run a supervised classification  on a landsat 8 (L8 OLI/TIRS COLLECTIONS LAND SURFACE REFLECTANCE ON-DEMAND) imagery. I created a signature file, then ran the Maximum Likelihood classification tool, but the output wasn't what I anticipated- which  was to get a classified raster. Instead I got a black and white raster with high to low values ranging  from 152-62.

I used the Image Classification Window (training samples, maximimum likihood, etc.) as seen:

I created the training samples with the draw polygon option, then used the Training Sample Manager to create he classes.  The classes were, Agriculture, Water, Urban, and Bare Earth. The values used were automatically generated by ArcMap, as seen below:

The parameters used were as follows:

The results were:

I tried running an unsupervised classification,  and that did s classification,  so I'm not exactly sure what I'm doing wrong.


texts: How to delete the existing polygon feature value inside the new polygon so that it borders the new polygon instead of overlapping?

texts:
Using arcGIS 10.7
Trying to create a polygon feature inside an existing polygon using the Create feature tool in Edit mode. How do I delete the existing polygon feature value inside the new polygon so that it borders the new polygon instead of overlapping? Is there a way to draw the new polygon so that it simultaneously creates the border and removes the existing polygon within?


texts: General terms for the Esri-terms "Layer" and "Feature Class"

texts:
I'm familiar with these terms used by ESRI:

Layer

Stores visual properties (color, symbols,…)
does not store data (attributes, coordinates)
links to a Feature Class (or raster data)

Feature Class

a spatial dataset (attributes, coordinates)
collections of points, lines, or polygons and their attributes
no symbols (color, …)
e.g. stored as shapefile on disk or feature class in a geodatabase

Exist for them terms that can be used for (nearly) all GIS (independent of the manufacturer)?

I think "layer" can be used for all GIS. But "Feature Class"?

I need the terms because I want to describe the functionality of a GIS independently of a manufacturer.


texts:
I work with several maps with this text:

and i try to change the second line only into "landUse" by using this code:

import arcpy
from arcpy import env

env.workspace = r"G:\desktop\Project"
for num, mxdname in enumerate((arcpy.ListFiles("*.mxd")),start = 1):
    print '\n',num,mxdname
    mxd = arcpy.mapping.MapDocument(
    r"G:\desktop\Project\\" + mxdname)
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        elm.text = elm.text.replace(u'aaa'[1], u'landUse')
    mxd.save()
del mxd

when i run it i get this unwanted result:

whereas i want this result:

i didn't find an answer in Changing part of text element using ArcPy?


texts: mapquest geocoding - i can not get the (lon,lat) returned from the API

texts:
I just started using the mapQuest Geocoding API and I created a code for searching addresses.
The following code alerts "OK" but I can not alert the (lat, lon) coords, even I tried with so many different ways.
Can you spot what is wrong? Does this means that the coords are not actually returned?
The alerts about the coords are not even alert "undifined"...

Thanks

the code :

function codeAddress(){
    var ad = document.getElementById("address");//get input from form
    $.ajax({
                type: "POST",
                url: "http://www.mapquestapi.com/geocoding/v1/address?key=Fmjtd******&outFormat=json&inFormat=json&json={location:"+ad+",options:{maxResults:1}}",                                              success: function(response) {
                    alert("OK");//alerts OK, but NONE of the followinng
                    alert(response.results[0].locations[0].location.street);
                    alert(response.results[0].locations[0].location.latLng.lat);
                    alert(response.results[0].latLng.lat);
                    alert(response.results[0].location.latLng.lat);
                    alert(responce.results[0].locations[0].latLng);

                    }})
}

//this is the html form
<input id="address" type="textbox">
<input type="button" value="Encode" onClick="codeAddress()">


texts:
I used the following line to try and add a vector shapefile to my qgis map, though it gives me "invalid layer" errors. Does anyone see something I don't?

iface.addVectorLayer(r"/Users/pgcseismolab/Desktop/jaden_pgc_archive/Wells/WD/Shapefiles/wd_14/wd14.shp","test","ogr")

I'm fairly new to python- do I have to define a variable before running this command? *on Mac OSX


texts: How to get index of specific location in IDL

texts:
I am trying to get the index of specific lat and lon in IDL.

I know how to do it in python.

tmp = tmp.sel(lat=lat_loc,lon=lon_loc,method='nearest')

But I can't find a way to do so in IDL.

Any suggestion or idea?


texts: Converting field to string with leading zeroes using Python Parser of ArcGIS Field Calculator?

texts:
I would like to create an expression using the field calculator in ArcMap to be able to convert an Int field to a string and add leading zeroes.

I am able to do it in a python window with an example

x=450
print str(x).zfill(4)
>>>0450

I have tried:

!Time!.str().zfill(4)

but that gives me invalid syntax errors.


texts: How to open overlay content from leaflet popup

texts:
I am currently working on an interactive leaflet map. I managed to display markers from a database via php (create an array). Now I want users to create markers with a click on the map. To add a marker with a click is not the problem, nor to delete the marker again. I have created two buttons in a popup (bind to the created marker). While one marker deletes the marker, the other button shall open a form.

I want this form to be displayed with jquery overlay function (in front of the map and the popup). I managed to open the overlay content, when the link is outside the leaflet-script. But I have no idea how to do open the overlay content from inside the leaflet-script. It just doesn't work!

Does anybody had kind of the same problem? I am not very sure if that will work at the end, I build my script step by step. But an overlay function would be very nice to have.

Hope this picture explains my problem.

<title>Marker Leaflet Map</title>
<meta charset="utf-8" />
<link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7/leaflet.css"/>

" alt=""/>


texts: OGR Geometries Error When Writing to Disk: 'Incompatible geometry for operation'

texts:
I am working on a piece of code that does a pretty complex computation using GPKG layer feature areas, I then calculate several new fields . Finally, when I try and write the results back to disk, I get an inconsistent (doesn't always appear) and nondescript error from what I assume is OGR: 'Incompatible geometry for operation'.
I am working with a LOT of features, so I am utilizing multi-threaded processing. All of my vector data is stored in zip files, so I use temporary directories to unzip, process and then re-zip the data. The line of code where I think I get the error looks like this:
dst_vector_disk_ds = ogr.GetDriverByName('GPKG').CopyDataSource(out_mem_ds, vector_file.replace('features', f'split_{job_index}'))
where vector_file is the original vector file name (string), out_mem_ds is the in-memory data set with all the processed features now neatly packaged away in their respective layer and job_index is the thread number referring to the multi-threaded processing.
I am using WSL (Windows 10), Anaconda3, Python 3.8.2 environment, GDAL version 3.0.2 and I have already used the flag ogr.UseExceptions().
Haven't been able to pinpoint the actual source of this error and since my Python code actually calls on a .c file module (compiled with Cython), de-bugging this has been very difficult.
The vector data is a layer of polygons. The threads write temporary outputs and then pass those back to the parent which merges all the resulting temporary output files to a permanent file.I have narrowed the source after a lot of manual de-bugging: The incompatible geometries come from the in-memory features that I am trying to write out to disk not having their original geometries associated. I have tried storing their geometries and re-associating them, but I ran into one of the classic ogr reference pitfalls/gotchas as mentioned in the GDAL Gotchas.


texts: Creating shp from txt files using ArcPy?

texts:
I have a txt file like this:

Point
0 347172.2 4269952.7 60.9090909091
1 337172.8 4278952.4 85.7334909091
2 344572.3 4547952.2 123.909673091
3 347532.8 8777952.8 456.467809091
END

The third column are some values i want to have in my points.shp. I wrote this script but it doesn't seem to work, meaning it doesn't pass the specific values of the third column

import arcgisscripting
gp = arcgisscripting.create()
# Set up inputs to tool
inTxt = arcpy.GetParameterAsText(0)
inSep = "."
strms = arcpy.GetParameterAsText(1)
# Run tool
gp.CreateFeaturesFromTextFile(inTxt, inSep, strms, "#")

What should i add also? I am using ArcGIS 10.2


texts: How to assign geometry to records in a table that contain an ID of a geometry point?

texts:
We are looking to map out 'incidents' at schools, but the incident table, which contains a few thousand records, only has a school number (plus the incident data). The table format is an Excel Spreadsheet with around 20,000 records.

We have schools (around 250 records/points) that share the same ID and of course the geometry for each school. The schools are a point SDE feature class, but can be converted to SHP or other geometry type if needed.

How can we assign the geometry from the schools to each record in the 'incidents' table?

I'd like to have something like a shapefile or other type of table. Keeping the excel file in its current format is not necessary, ie. we're converting the excel file to a feature class / shapefile / etc.

Any software option is on the table.


texts: Reduce regions on projected (equal-area) data in Google Earth Engine

texts:
I am working with a global raster of soil carbon and a feature collection of ecosystems. I would like to average the carbon in each ecosystem using ReduceRegions in Google Earth Engine, but I need to run this calculation in an equal area projection like EckertIV or Mollweide.
My code works  fine when I don't reproject the data, but the values are incorrect at a global scale. I have been trying to reproject both the raster data and transform the feature collection to Mollweide projection and then use Reduce Regions (as seen in the code below), but I continue to get a "User Memory Capacity Exceeded" error when exporting the results as a csv which did not happen prior to reprojecting.
var wkt = ' \
  PROJCS["World_Mollweide", \
    GEOGCS["GCS_WGS_1984", \
      DATUM["WGS_1984", \
        SPHEROID["WGS_1984",6378137,298.257223563]], \
      PRIMEM["Greenwich",0], \
      UNIT["Degree",0.017453292519943295]], \
    PROJECTION["Mollweide"], \
    PARAMETER["False_Easting",0], \
    PARAMETER["False_Northing",0], \
    PARAMETER["Central_Meridian",0], \
    UNIT["Meter",1], \
    AUTHORITY["EPSG","54009"]]';

//Define the Mollweide Projection
var Proj_Mollweide = ee.Projection(wkt);

//Reproject SOC Raster to Mollweide at 250m resolution
var SOC_1m_tCHa_Mollweide = SOC_1m_tCHa.reproject({
  crs: Proj_Mollweide,
  scale: 250
});

// Transform Biomes into Mollweide
var transform_mollweide = function(a) {
  var transformed_mollweide = a.transform(Proj_Mollweide, 1);
  return transformed_mollweide;
};
var Biomes_Mollweide = Biomes.map(transform_mollweide);

//Define the reducers
var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.count(),
  sharedInputs: true
});

//Reduce the SOC data within each biome polygon
var SOC_by_Biome = SOC_1m_tCHa_Mollweide.reduceRegions({
  collection: Biomes_Mollweide,
  reducer: reducers
});

Here is the link to the code I have been working with: https://code.earthengine.google.com/3db50c3c74fe41428c62c3f7193e131f
Is there a way to complete this calculation with the reprojected data in Google Earth Engine?


texts: QGIS Print Layout Vertical Graduated or Continous Legend Bar

texts:
I would like to replicate a legend bar design from a map that I assume comes out of ArcGIS:

Is this possible in QGIS without using external software like illustrator?


texts: Searching open source Terrain and Satellite maps for Leaflet

texts:
For our project, we are exploring Terrain and Satellite maps using Leaflet.js

While exploring the Internet got this OpenTopoMap. However, not sure whether this is a Terrain Map.

For Satellite got this Esri.WorldImagery. However, as the map source is https://server.arcgisonline.com don't know whether these are open source maps or paid maps.

Also, is there any open source basemap like OSM for Terrain and Satellite maps?


texts: (Uncaught SyntaxError: Unexpected token < in JSON at position using openlayers js)

texts:
I am trying to insert a map with WFS request from GeoServer with OpenLayers tool but I get this error:

Uncaught SyntaxError: Unexpected token < in JSON at position 0
    at JSON.parse (<anonymous>)
    at Jl (ol.js:342)
    at wm.k.Pa (ol.js:342)
    at T.<anonymous> (ol.js:340)

and this is my code:

    <link rel="stylesheet" href="https://openlayers.org/en/v4.1.1/css/ol.css" type="text/css">

I don't know where is the problem because I didn't found anything with JSON problems and OpenLayers.


texts: Control the colours of shapefiles loaded to Geoserver?

texts:
I have a set of shapefiles loaded to GeoServer, which I later deliver as a WMS. All the shape files are rendered in the same neutral grey, which makes it quite difficult to tell apart which shape file I'm visualizing. Can I assign colours of textures to each polygon?
A very similar question was posted earlier. But, since the user relied on ESRI technology, all answers focused on ESRI solutions; not GeoServer.
Control the colors of a shapefile
The question linked to mark mine as duplicate is a different question; even if both answers happen to be in the same manual. I want to give a colour to each shapefile. While the other question was about assigning colours based on the value of a field.


texts: Exporting mxd using ArcPy/Python without opening it?

texts:
I would like to be able to update a load of pdf and jgeg maps once the data has been updated. Currently people have to open each mxd and export the map in each format. I've written a short python script that will export the current mxd. This works well:

##export current mxd to pdf and jpeg
mxd = arcpy.mapping.MapDocument("CURRENT")
mxd.save()
filename = mxd.filePath.split("\\")[-1].split(".")[0]
path = '\\\\main.corp.local\\E-GB$\\Home\\4\\J042\\Documents\\PWISE' + "\\" +filename
arcpy.mapping.ExportToPDF(mxd, path + '.pdf')
arcpy.mapping.ExportToJPEG(mxd, path + '.jpg',"PAGE_LAYOUT", 400,500,300)

This doesnt save all that much time though because someone still needs to open all the mxds and run it. I've therefore written this to run on each mxd:

##export current mxd to pdf and jpeg
import arcpy, os
#mxd = arcpy.mapping.MapDocument("CURRENT")
for root, dirs, files in os.walk(r"\\main.glb.corp.local\EP-GB$\Home\ABZ\4\J0422414\Desktop\EXPORTMXD"):

        for f in files:
            if f.endswith(".mxd"):
                mxdfile = root + '\\' + f
                mxd = arcpy.mapping.MapDocument(mxdfile)
                mxd.save()
                filename = mxd.filePath.split("\\")[-1].split(".")[0]
                path = '\\\\main.corp.local\\E-GB$\\Home\\4\\J042\\Documents\\PWISE' + "\\" +filename
                arcpy.mapping.ExportToPDF(mxd, path + '.pdf')
                arcpy.mapping.ExportToJPEG(mxd, path + '.jpg',"PAGE_LAYOUT", 400,500,300)
            else:
                print f

The weird problem here is that when there is a table in the layout view of the lap (as in attribute table  menu > Add Table to Layout) it will not show - there is just a white space. This problem occurs both on the pdf and the jpeg

Does anyone know either if there is a solution to this export bug, or if I can adjust my script to open and close each mxd, running the relevant part of the script from within it? because it seems to work if the mxd is open


texts:
I need to apply Polygonize tool to several hundred shapefiles with Graphic Modeler tool. However, the Polygonize tool empties the attribute table - and I need to retain at least one of the attributes.

Is it possible to do in a Modeler?


texts:
In arcgis office you have a functionality to add a map to an excel.

Now, i want to give an option to export an excel from my application(javascript, .net, or anything else) to allow the user to download an excel with the map already embedded.

It is possible? How?


texts: Transform bookmarks to data driven pages in ArcGIS Desktop?

texts:
We have this large project in ArcGIS with 30 maps with different extents using Bookmarks to pan between these maps/pages.

Now we would like to use Data driven pages instead. The question is, is it possible to convert bookmarks to data driven pages or do we have to create a lot of polygons covering each map extent?


texts: Polygon label parallel to line but inside polygon

texts:
I am trying to add labels to a polygon layer hosted on GeoServer with CSS styling, and I would like the labels to be on the boundary running parallel, but always inside of the polygon.
For instance if there is a horizontal boundary, I want the label for the top polygon to be above the line, and the label for the bottom polygon to be below the line.
I am using the CSS plugin for GeoServer, and using this styling
[@sd >=0] [@sd <= 1000000] {
  font-family: "DejaVu Sans Book";
  font-size: 18;
  font-fill: #C3282D;
  label: [DIST_NM];
  label-anchor: 0.5 0.5;
  label-offset: -10;
  label-max-displacement: 50;
  halo-radius: 1;
  halo-color: white;
}

but the label offset of -10 appears to be relative to the line (above or below) and not the polygon (inside or outside).
In the scenario described above with this styling, the labels for the top polygon and the bottom polygon are both drawn below the line, so it is unclear what polygon the label belongs to.
For additional context, I was able to get this kind of styling with mapnik using this kind of rule
<Rule>
    <MaxScaleDenominator>1000000</MaxScaleDenominator>
    <MinScaleDenominator>0</MinScaleDenominator>
    <TextSymbolizer dy="-5" size="18" face-name="DejaVu Sans Book" spacing="500" placement="line" vertical-alignment="top" fill="#C3282D" halo-fill="white" halo-radius="1" >[DIST_NM]</TextSymbolizer>
</Rule>


texts:
I am practicing QGIS python API and I made standalone script to create .qgs file and import some layers from postgis database. Now I need to zoom to imported layer but I saw that it is not possible via iface.setActiveLayer(vlayer) and  then iface.zoomToActiveLayer() because I cannot use iface for standalone script.

Is there any other way to zoom to imported layer?

This is what I have done so far:

import subprocess
subprocess.call([r'batch.bat'])
import qgis.core
print('Importovan')
import glob
import os
import sys
from qgis.core import (QgsProject, QgsRasterLayer,QgsVectorLayer,QgsApplication,QgsCoordinateReferenceSystem)
from qgis.gui import QgisInterface,QgsMapCanvas
from PyQt5.QtCore import QFileInfo
from PyQt5.QtSql import *
from qgis.core import QgsDataSourceUri
import qgis.utils
strProjectName = "my_project.qgs"

# start the qgis application
QgsApplication.setPrefixPath(r"C:\OSGeo4W64\apps\qgis", True)
qgs = QgsApplication([], False)
qgs.initQgis()

# start a project
project = QgsProject.instance()

selectedcrs="EPSG:4326"
target_crs = QgsCoordinateReferenceSystem()
target_crs.createFromUserInput(selectedcrs)
project.setCrs(target_crs)

canvas=QgsMapCanvas()
iface=QgisInterface.QgsMapCanvas()
project_path = strProjectName
###############################

db = QSqlDatabase.addDatabase('QPSQL')
# check to see if it is valid
if db.isValid():
    print ("QPSQL db is valid")
    # set the parameters needed for the connection
    db.setHostName('localhost')
    db.setDatabaseName('request')
    db.setPort(int(5433))
    db.setUserName('postgres')
    db.setPassword('postgres')
    uri = QgsDataSourceUri()
    uri.setConnection("localhost", "5433", "request", "postgres", "postgres")
    #open (create) the connection
    if db.open():
        print ("Opened %s" % uri.uri())
        uri.setDataSource ("rodno_razvrstani_podaci", 'statistika_2017', 'wkb_geometry')
        vlayer=QgsVectorLayer (uri .uri(False), 'statistika_2017', "postgres")
        project.addMapLayer(vlayer)

    else:
            err = db.lastError()
            print (err.driverText())
db.close()

#This is where script breaks
qgis.utils.iface.setActiveLayer(vlayer)
#iface.zoomToActiveLayer()

###############################

# write the project file
project.write(project_path)

# stop qgis
qgs.exitQgis()


texts: Generate in oracle sql developer a point with specific SRID from two variables

texts:
I have created quite often (in posgis) points using two variables from the table (LAT,LON) and setting the desired SRID.
I am trying to do something similar now with Oracle SQL developer but I keep getting errors when I try to display this as a layer in geoserver.
This is what I try:

SELECT
  ID,
  CO_NAME,
  ADJUSTMENT_FACTOR,
  MDSYS.SDO_GEOMETRY(null,2100,SDO_POINT_TYPE(24.004632439599,24.004632439599,NULL),NULL,NULL) AS GEOM
 FROM CENTRAL_OFFICES

Instead of the variables LAT and LNG I have just put two random numbers. Moreover I want to ask how is possible to view in a readable format the generated geometry. Something like: ST_ASTEXT of postgis.

EDIT

I have written this query:

SELECT ID,
  CO_NAME,
  ADJUSTMENT_FACTOR,
  PORT_LOCATION,
  VDSL_ENABLED,
  SDO_UTIL.TO_WKTGEOMETRY(SDO_CS.TRANSFORM(
    SDO_GEOMETRY(2100,4326,
      SDO_POINT_TYPE(CAST(REPLACE(CENTRAL_OFFICES.LAT, '.', ',') AS FLOAT), CAST(REPLACE(CENTRAL_OFFICES.LNG, '.', ',')
        AS FLOAT),NULL),NULL,NULL),2100))    AS  geom
  FROM CENTRAL_OFFICES

This returns a table with a point geometry in the desired SRID (2100).


texts: I need to generate routes that will be performed in long vehicles. The path can not have curves 90 degrees because the streets are narrow.


texts: I need greenspace data for Brighton and Hove and Dudley areas. Does anyone know where I can get this? I have tried to extract it from OpenStreetMaps but the areas are too large and it crashes before the data is downloaded.


texts: Why does Georef'd image show as a large black area

texts:
Using QGIS2.4.0 I just added 2 maps to an existing project. Both are scanned images of pre-1800 plans. One is a jpg file, I georef'd and added to QGIS no problem and plan displays in the project.
I carried out exactly the same process with the other plan, a tif file, and the process completed  but the plan just displays as a black rectangle, can anyone advise a very novice novice what might have gone wrong please.

Could it be because the jpg was a colour image, but the tif is a greyscale? I've had problems with greyscale images before.


texts:
I have a vector layer with millions of polygons making continuous coverage. I need to classify them according to their shape. I am already using several shape indexes from landscape ecology like compactness (4piA/P^2), mean width (2A/P), shape number (P/sqrt(A)), I also saw this answer to Calculating roundness/compactness of polygon?

My problem is that all these metrics are using some ratio of area and perimeter only. Even the Fractal Dimension index is using only area and perimeter (2ln(0.25P)/ln(A)). But how can I distinguish two polygons with the same area and perimeter but absolutely different shape? Like this branched polygon A:

which I tried to draw with the same area and perimeter as curved strip B. All my known indexes will be the same for them. But for me it's very important to differentiate simple strips (including curved like new moon) from complex branched shapes.

I intentionally show the polygon B as a curved strip and not a straight strip because I am aware of Related Circumcircle index which detects straight elongated shapes but my polygons may have also the same circumcircles. Even if I construct Convex hull and calculate a ratio of areas Apolygon/Aconvex, it may be very similar here.

So, how can I clearly distinguish branched polygon A from polygon B in vector data automatically? (Converting them to raster would require extremely small cell size, enormous dataset and a lack of memory, so it is not possible). Are there other shape indexes which include other parameters? Ideally, the method would distinguish not only clearly branched polygons but even C and D:

My only idea is to construct the convex hull then erase the polygon from it's convex hull and count the number of (big) pieces it leaves (erasing polygon by polygon and not the entire layer). This could show the border complexity.

I welcome mathematical solutions/algorithms, which I would later implement in Python.


texts: How to create a dropdown menu from another dropdown menu for the ArcGIS Collector App?

texts:
I'm creating a point Feature Class in ArcMAP for the Collector App and I'm trying to create an option to select an attribute, and based on that selection another set of attributes would be available to choose from. I believe it's called a cascading dropdown list.
The Field is "Identity" and from there I set up a Domain where the options are "Utility, Drainage Features, Geology, etc...". For example, when they choose "Utility" the set of options should be "Electrical, Gas, Water, Telecom". I tried to use subtypes but then it just changed the dropdown menu to the Utility options and there was no longer access to the "Utility, Drainage Features, Geology, etc..." attributes.
Is there a way to do this?
Otherwise I may have to just create another Field called "Type" with all of the specific possible attributes.


texts:
Has anyone tried linking a Sharepoint Excel spreadsheet with CartoDB?

My spreadsheet has a URL like it would from Google Drive or Dropbox and it is shared with editing and viewing privileges for everyone (public),  however it keeps erroring when I try to import.

error:

Error 1001

Download error

What about...

The remote URL returned an error. Please verify your file is available at that URL.

Any ideas? I would like to use Excel as the program for my team to edit my data and I want it synced to CartoDB.


texts: Plotting latitude on (Mercator projection) SVG representing country map

texts:
I want to convert (google sourced) longitude and latitude coordinates to x and y to place points on a SVG map
I can get the longitude converted to X but I can't nail latitude.
http://jsfiddle.net/chrisloughnane/an3BZ/17/

red dots = place holders
green dot = calculated position from longitude/latitude

I have read so much about Mercator projections and other projection system API's I am now completely confused. I also followed Proj4JS library threads but couldn't find an example that was close to my task.
I attempted to emulate @Raphael Wichmann and @betamax solution, unfortunately I got something wrong.
Could someone have a look at my jsfiddle and see if it's an obvious mistake?
The second test function secondconvert(latitude, longitude) (bottom of javascript pane) provides me with a reasonably accurate x coordinate i.e. if I manually enter the y coordinate b.ylat change to 265 it covers my left red dot place holder nearly perfectly.
I think I'm close.
Note: the longitude/latitude bounds I used for my SVG map are best guesstimates, worked surprisingly well for the longitude calculation.

Original SVG from here.


texts: Cutting a hole out of a LAS file to omit that area?

texts:
I have a LAS dataset of a forest stand and a chunk of it has been recently clearcut.  I would like to update the LAS to reflect this so I can model the area but showing the clearcut hole in the canopy.

Is there a way to either remove those LAS points to leave a void in the middle of the dataset, or alternately I guess I could reclassify the points into one of the unassigned classes and then ignore that class of points?

I would like to accomplish this using QGIS and would like to alter the dataset or the underlying LAS files, whichever method will accomplish removing the cutblock from the display of the LAS points.


texts: Making PostGIS polygons ready to use in ArcGIS

texts:
I've got a large dataset in PostGIS which I've made available to an indirect colleague who's using ArcGIS. However, he can't load most of the features which have the polygon geometry type.
I've already tried using ST_MakeValid and furthermore, I've validated all features with ST_IsValid which they are.
The error message he's getting in ArcGIS is a "Shape integrity error". Does anyone have an idea what it could be? Furthermore, does anyone know a good workflow in which I can troubleshoot this myself? Since I don't have the ArcGIS software I can't verify my solutions.


texts: The example provided in http://api.geoext.org/1.0/examples/tree.html shows all overlayes are checked but how do I make these overlayes unchecked while page is loaded?


texts:
I'm working with rasters, there are 3 folder (input, query, result). In the input folder are the basic rasters. The user can choose between these rasters, the code can collect the rasters into the query folder, and after it will make a mean raster into the result folder. But after the second run i need to clear the query folder, because it stores the first run rasters and it gives a wrong result. How can i delete these rasters from the query folder, after every run?

 import os, sys
import shutil
import arcpy
from arcpy import env
arcpy.CheckOutExtension("spatial")
from arcpy.sa import *
arcpy.env.overwriteOutput = True
start_date = arcpy.GetParameterAsText(0)
end_date = arcpy.GetParameterAsText(1)
env.workspace = "c:/temperature/input"
rasters = arcpy.ListRasters("*","GRID")
for raster in rasters:
    if int(raster) >= int(start_date) and int(raster) <= int(end_date):
        arcpy.CopyRaster_management((raster), "c:/temperature/query/" + (raster) + ".tif", "", "", "", "NONE", "NONE", "", "NONE", "NONE")
env.workspace = "c:/temperature/query"
eredmeny = "c:/temperature/result/mean.tif"
rasterses = arcpy.ListRasters("*","TIF")
result = arcpy.gp.CellStatistics_sa((rasterses), result, "MEAN", "DATA")
exit()


texts:
I was wondering if there is a way to load a local tile cache (ie: MB Tiles etc.) of imagery that is sitting locally into QGIS?  I am hoping to use some tablets for field verification and would like to have ortho imagery to assist with the verification.  Where I will be doing the field work I will not have wifi or cell coverage, so I will need local connections.  I am hoping to use tiles over source imagery to minimize the large storage requirements needed for the orthos, and to increase the performance on a tablet.

Thanks...


texts: Rasterise 3D polyline using z-values to interpolate cell values

texts:
Using ArcMap 10.5.1.
I have a ZM polyline composed of 3x vertices, with the vertex z-values representing water level at wells.  I am looking for a tool to rasterize this polyline, using the z-values to linearly interpolate the raster cell values between the wells (the ultimate aim is to then use these values with the Euclidian Allocation tool and other inputs to calculate the water table level).

I am aware that such a facility exists in GDAL (see: Rasterise heights from a 3D linestring) but I cannot find an ArcGIS tool to do the same. The examples I can find (e.g. Interpolate Shape) are for adding z-values to a 2D polyline from an existing raster.


texts: QGIS "Execute SQL" don't accept more than 1 input table

texts:
I have some problems using the "Execute SQL" algorithm (in General Vector Tool) in the Processing Toolbox.
No problem when I've got only 1 input table, but when I try to select 2 tables as input data I get an error: virtual: Referenced table input2 in query not found!

For background information my SQL expression (works fine with add virtual layers):

select makepoint((st_line_locate_point(p.geometry, a.geometry)*st_length(a.geometry)) /*:real*/,"Z" /*:real*/,2154) as geometry
from 'input1' as a, 'input2' as p


texts: get center coordinates of a buiding(way)

texts:
Using overpass-turbo I got data like this:

{
  "type": "way",
  "id": 292833530,
  "nodes": [
    2963934681,
    924717449,
    2963934682,
    2963934683,
    2963934684,
    2963934685,
    2963934686,
    2963934681
  ],
  "tags": {
    "addr:city": "Ivry-sur-Seine",
    "addr:housenumber": "134",
    "addr:postcode": "94200",
    "addr:street": "Avenue Danielle Casanova",
    "amenity": "post_office",
    "atm": "yes",
    "building": "yes",
    "change_machine": "yes",
    "moneo:loading": "yes",
    "name": "Ivry principal",
    "operator": "La Poste",
    "phone": "3631",
    "ref:FR:LaPoste": "17313A",
    "source": "data.gouv.fr:LaPoste - 01/2013",
    "stamping_machine": "yes",
    "wheelchair": "yes"
  }
}

Is there a way to get approxmatly the center of this building?
I use

  http://api.openstreetmap.org/api/0.6/node/ + one of the node

to get coordinate for un point but I would like something more in the center of the building.


texts: Finding Cubesat Data for remote sensing dissertation?

texts: Where can I find CubeSat data so that I can compare it to other satellite data for my dissertation?


texts: How to show marker information into a tab

texts:
I have a lot of information that I want to show on a Leaflet Popup. For these type of markers (I have four different layers and each marker show different information) I have three elements:

General description
Images
Video

But I want to show them as tabs inside the Popup. Something like these:

To add each icon I use these method:

    L.marker([20.683, -88.568], {icon: tiloIcon}).bindPopup('
        Chichén Itzá
Cultura
Maya
 Descripción
Es uno de los sitios arqueológicos más importantes de la Península de Yucatán. Es la ciudad más impresionante del periodo Clásico maya y fue inscrita en la lista del Patrimonio de la Humanidad de la UNESCO.
Ver ficha completa
            ').addTo(sitios),

Additionally, I want that the title tab disappear and show it at the top of the popup

Do you know any way to put the information on tabs?


texts: Finding '_' without treating it as a wild character

texts:
Does anyone know how to select records containing an underscore without the query treating it as a wild character?

For instance, I'm trying to select records in a text field ending with '_' to then right trim the underscore off the record.


texts:
I am trying to run a MaxEnt in ENMeval to select the best performing model with differing regularization multipliers and feature combinations.
I am not that strong in RStudio but have been advised that this is the best approach I have been following an online tutorial but have come across an error in my code:
sp_df <- read.csv(file = "Season1.csv", header = T);

sp_df <- cbind(sp_df, rep.int(1, length(nrow(sp_df)))); #Adds another column indicating these are presence points

colnames(sp_df) <- c("Species", "X", "Y", "Response");

#Get environmental variables

preds <- list.files(pattern = ".asc");

preds <- stack(preds);

eval = ENMevaluate(occ=sp_df[,c('X', 'Y')], env = preds, method= 'randomkfold',

kfolds=10, parallel=TRUE, numCores = 12, fc=c("L", "Q", "LQ"), RMvalues=seq(0.5, 2, 0.5), rasterPreds=T)

*** Running ENMevaluate using maxnet v.0.1.2
*** Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘nlayers’ for signature ‘"character"’

How can I fix this error? I'm confused as I don't understand the language in this error or what it is referring to?


texts: PostGIS/QGIS Error geometry has Z dimension when snapping

texts:
I am working with QGIS 2.14 and PostGIS 2.3 (for a PostgreSQL 9.5), and I am hitting a strange error.

I have a SHP with some areas, which I would like to include in a PostGIS table (MULTIPOLYGON), manually because I want to add some attribute informations. It works fine when I simply create a polygon, but I want to follow the SHP's polygons, so I activated the snapping option, but when trying to save, I get the error Geometry has Z dimension but column does not.

A logical explanation for me would be that the snapping tool create a Z attribute (even if it has a NULL or 0 value), and the error would therefore occur when trying to save this coordinates in a 2D geometry... But that's just a theory.

Does anyone know what is happening when using the snapping ? Is there a consideration of a Z dimension ? And if so, would it be possible to deactivate it ?


texts: Convert CAD based on known Local Origin Point

texts:
I have CAD files which were created in a "Staten Island local datum". There is no spatial information associated with these files, just annotation mentioning the datum. After researching this datum, I discovered that it's origin is located at N 137,190.02, E 1,967,746.81, according to the NAD_1927_StatePlane_New_York_Long_Island_FIPS_3104 coordinate system.

Given this info, how do I get these files to line up with the rest of my data? Do I need to create a custom projection?


texts: What do horizontal (H) and vertical (V) polarization mean for satellite radar?

texts: I'm reading about synthetic aperture radar (SAR), and one of the basic concepts is that the transmitted and received signals can have "horizontal" (H) or "vertical" (V) polarization, but I can't find anywhere what that means in terms of the real-life physical direction.  If the radar were on top of a building and transmitting a wave parallel to the surface of the earth, then I can understand that horizontal and vertical would mean relative to the surface of the earth.  But if the radar is on a satellite, then the direction of propagation of the wave is more like perpendicular to the surface of the earth, so both of the polarizations would be nearly horizontal in that sense (and in the case where the radar is looking directly downward (zero "incidence angle"), both polarizations would be exactly horizontal relative to the surface of the earth).  So what do horizontal and vertical refer to in this context?  Is vertical the direction that's in the plane containing the direction of the satellite's movement and the direction of the wave (and is perpendicular to the direction of the wave), and horizontal the direction that's perpendicular to that plane?


texts: Using the OpenStreetMap Service in PyQGIS

texts:
i am trying to learn how to use OSM Service in QGIS using Python.

I wrote this script:

import os
from osgeo import gdal
from qgis.core import QgsRasterLayer, QgsDataSourceUri

project = QgsProject.instance()
project.clear()

xml = """
<GDAL_WMS>
    <Service name=”TMS”>
        <ServerUrl>http://tile.openstreetmap.org/${z}/${x}/${y}.png</ServerUrl>
    </Service>

    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>19</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>

    <Projection>EPSG:3857</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache />
</GDAL_WMS>"""

vfn = "/vsimem/osm.xml"
gdal.FileFromMemBuffer(vfn, xml)

rlayer = QgsRasterLayer(vfn, "OSM")
if not rlayer.isValid():
    print("Raster failed to load!")
else:
    project.addMapLayer(rlayer)

But when I run It, It print "Raster failed to load!".
Where did I make the mistake?


texts: Default or wrong behavior with multiple clicks in intersection zone in OpenLayers?

texts:
We use Vue.js and OpenLayers 4.6.5 in our project. I've noticed the following behavior - when I have two intersected polygons with some styles and click first time inside the intersection area (look at the arrow pointer), the first polygon is selected. When I click second time in the same area, the other polygon is selected.

I've checked the official documentation and found nothing about this. Is this default behavior? I've also checked my code and I can't see that I did something that could produce it.

My code is as following:

markObject (mark) {
  if (!mark) {
    if (this.lastSelectedFeature) {
      const style = this.lastSelectedFeature.get('originalStyle')
      this.lastSelectedFeature.setStyle(style)
    }
    this.map.un('select', this.onMarkObject)
    if (this.markSelection) {
      this.markSelection.getFeatures().remove(this.lastSelectedFeature)
      this.map.removeInteraction(this.markSelection)
    }
    return
  }

  if (!this.markSelection) {
    this.markSelection = new Select({
      condition: condition.click,
      layers: [this.vectorLayer]
    })
    this.markSelection.on('select', this.onMarkObject)
  }

  this.map.addInteraction(this.markSelection)
},
onMarkObject (event) {
  if (event.selected.length === 0) {
    if (this.lastSelectedFeature) {
      const style = this.lastSelectedFeature.get('originalStyle')
      this.lastSelectedFeature.setStyle(style)
      return
    }
  } else {
    const clickCoordinates = event.mapBrowserEvent.coordinate
    const allFeatures = this.vectorLayer.getSource().getFeatures()
    const drawnFeatures = allFeatures.filter(x => x.get('graphics'))

    for (const feature of drawnFeatures) {
      if (feature.id_ !== event.selected[0].id_) {
        const isIntersects = feature.getGeometry().intersectsCoordinate(clickCoordinates)
        if (isIntersects) {
          console.log('yes')
        } else {
          console.log('no')
        }
      }
    }

    if (event.selected[0] === this.lastSelectedFeature) {
      console.log('click the same')
    } else {
      console.log('click another')
    }
  }

  if (event.selected && event.selected[0].get('graphics') && event.selected.length > 0) {
    if (this.lastSelectedFeature) {
      const style = this.lastSelectedFeature.get('originalStyle')
      this.lastSelectedFeature.setStyle(style)
    }

    this.lastSelectedFeature = event.selected[0]
    // Highlight currently selected feature
    if (event.selected[0].getGeometry().getType() === 'Point') {
      event.selected[0].setStyle(markPointStyle)
    } else {
      event.selected[0].setStyle(markStyle)
    }
  }
},

The screen shots:

How can I disable this behavior and write custom code?


texts: Condition with using geometry OLD and NEW

texts:
Is there a posibility to use a condition with OLD and NEW parameters?
Using postgres 12, and pgAdmin
Geometry in Qgis
Trigger
CREATE TRIGGER abc_trigger
AFTER INSERT OR UPDATE OR DELETE
ON "Table1" FOR EACH ROW
EXECUTE PROCEDURE abc();

Funcion
CREATE OR REPLACE FUNCTION abc()
RETURNS trigger AS
$BODY$
DECLARE
BEGIN
CASE TG_OP
     WHEN 'INSERT' THEN
     WHEN 'UPDATE' THEN
           IF OLD.stat <> NEW.stat OR OLD.geom <> NEW.geom THEN
                NEW.sum = (
                    SELECT COUNT(ot.id)
                    FROM other_table AS ot
                    WHERE ST_Intersects(st.geom, NEW.geom)
            );
         IF OLD.stat <> NEW.stat THEN
              UPDATE Table SET "sumary" = NEW.sumary WHERE Table."stat"=NEW.stat ;
         END IF;

    WHEN 'DELETE' THEN
END CASE;
RETURN NULL;
END;
$BODY$
  LANGUAGE plpgsql

What I need to do.
I`m count value for all table when table is UPDATE. But I want to do it only when row change stat or change geometry...
Does anybody know how to do this?


texts: Export spatial overlay from R package PBSmapping

texts:
After realising that the GEOS-library in Grass and the R package rgeos is incredibly slow, i found the packe "PBSmapping" to do the same thing, a Union of two Polygons and to cut the parts with no overlay-areas.

The problem is to export the created Polyset within R.
I use the following code:

#require
require("PBSmapping")
require("maptools")
require("rgdal")

# import .shp
area <- readShapeSpatial("area")
area.ps <- combinePolys(SpatialPolygons2PolySet(area))

area_WAL <- readShapeSpatial("area_WAL")
area_WAL.ps <- combinePolys(SpatialPolygons2PolySet(area_WAL))

# UNION
union.area <- combinePolys(joinPolys(area_WAL.ps, area.ps,"UNION", maxVert = 1.0e+06))

#export polyset as shape
export <- PolySet2SpatialPolygons(union.area)

... and get this error message:

In PolySet2SpatialPolygons(union.area) :
unknown coordinate reference system

As CRS I've used UTM and GK, neither of both worked.

Is there another way to export a Polyset-Class?


texts: Is there a way to get embedded output in the response using the qgis wps client.  It seems to me like it automatically asks for reference by url.


texts:
I'm trying to make a postgis table ready for pgrouting. When I run :

SELECT pgr_createTopology('network.VEJM_BR_SL', 0.001, 'geom', 'id');

I get an error:

  "Geometry has Z dimension but column does not"

It seems related to the "geom" column being the type "LinestringZ". I believe it means that it is with Z values. I've tried to use ST_Force2D like this:

SELECT pgr_createTopology('network.VEJM_BR_SL', 0.001, ST_Force2D('geom'), 'id');

but it does not work. Any suggestions?


texts:
I have a GeoTiff file with WGS84 datum projection and I'm trying to reproject it using gdalwarp. Using this command

gdalwarp -t_srs '+proj=tpers +azi=19 +tilt=30 +lon_0=-109.5 +lat_0=23 +h=1274199.4 +wktext' atnight.tif atnighttilted.tif

returns this

ERROR 1: Too many points (10000 out of 10000) failed to transform,
unable to compute output bounds.

I'm using GDAL version 1.10.0. This is the contents of the .tfw file:

0.0066667901
0.0000000000
0.0000000000
-0.0066669136
-32.8672753199
72.8560317049

I've just started using gdalwarp.
Thanks in advance,
Tom

The output of atnight.tif using gdalinfo is:

Driver: GTiff/GeoTIFF
Files: atnight.tif
       atnight.tfw
Size is 9123, 4957
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-32.870608714976186,72.859365161672642)
Pixel Size = (0.006666790125743,-0.006666913589395)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2012:11:26 15:13:57
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS6 (Macintosh)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -32.8706087,  72.8593652) ( 32d52'14.19"W, 72d51'33.71"N)
Lower Left  ( -32.8706087,  39.8114745) ( 32d52'14.19"W, 39d48'41.31"N)
Upper Right (  27.9505176,  72.8593652) ( 27d57' 1.86"E, 72d51'33.71"N)
Lower Right (  27.9505176,  39.8114745) ( 27d57' 1.86"E, 39d48'41.31"N)
Center      (  -2.4600456,  56.3354198) (  2d27'36.16"W, 56d20' 7.51"N)
Band 1 Block=9123x115 Type=Byte, ColorInterp=Red
Band 2 Block=9123x115 Type=Byte, ColorInterp=Green
Band 3 Block=9123x115 Type=Byte, ColorInterp=Blue


texts:
I want to do cubic spline interpolation with QGIS (3.6), but:

Why the result of interpolation is not a polygon?

And interpolation don't cover all country?

(see attached picture)?


texts: Converting Python Script back to a Graphic Model in QGIS Graphic Modeler

texts: I created a model in the QGIS Graphic Modeler and exported it as a python script. I made a few changes in the python script and now want to export it back to the model. Is it possible?


texts: How can I import pyqgis API into standalone Python script?

texts:
I want to create .qgis file from standalone Python script and I need to import pyqgis (qgis.core ). I tried commands:

set PYTHONPATH=C:\Program Files\QGIS 3.4\apps\qgis-ltr\python and set PATH=C:\Program Files\QGIS 3.4\bin;C:\Program Files\QGIS 3.4\apps\qgis-ltr\bin;%PATH%

and then in a Python console I tried

import qgis.core

and then I get error saying:

ModuleNotFoundError: No module named 'PyQt5.QtCore'


texts:
I want to write software that merges groups of (mainly) adjacent polygons. For example, merge the US counties into sales regions. Then, display either the boundary of these sales regions or fill them with some color indicating, let's say, sales numbers.

Preprocessing and merging the given shapes in a GIS editor is not an option, the sales groups are not given in advance.

The GIS component I use has a function for merging polygons, but it is awfully slow. Reducing the number of vertices didn't help much. Neither having the common segments of the polygons use the exact same vertices.
But given the fact that the line segments of the polygons are duplicated in a shapefile (which I use at the moment as data source), I can imagine that the algorithm might have a lot of work to do.

Are there algorithms and/or data structures that can handle this problem more efficiently?

I thought of maybe using Geo Markup Language or TopoJSON, which both have the ability to define segments which can be referenced to make up polygons.

Using that, is there software which optimizes existing polygons in such way that common boundaries are factored out into one representation and references to it?


texts: Is there a way to change the selection in the Editor Attributes window, while maintaining current selection?

texts:
I was wondering if anyone knew of a way to alter the selection within the Editor Attributes window while maintaining the current selection set.  In the image below, I've manually selected an item in the selection set so that attributes of the item are visible.

I was hoping to programmatically alter the selected item so that attributes of the newly selected item it would be visible in the attributes window.  I've seen the AttributeWindow Interface and the TableSelection Interface but these don't seem to do what I'd like, unless I've missed something. Is there another interface that I'm missing?


texts: ArcObjects(Java) - register table to geodatabase ESRI 10.3

texts: We are in need of connecting to the Oracle SDE schema for registering Tables to the sde schema using only ArcObjects(Java) is it possible to do that with no dependency. Assume the geodatabase was created from ArcDesktop tool installed on the WIN OS. However our need is to use ArcObjects Java SDK to connect to the sde schema from Linux OS.


texts: Read file form txt and perform operation by name in geopandas

texts:
I have a txt that contains in each line a sentence like this:

INT_ASTIK_EAS
ER_ASTOM_ASTIK
and so on.

I want to ask for a directory of the main folder that contains the subfolders that each contain the ASTIK.shp, then the EAS.shp, then the ASTOM.shp, etc

INT means to perform an intersection with this code where in this case should involve ASTIK and EAS shapefiles.

inte_s=gpd.overlay(ASTIK,EAS,how='intersection')

and check if it has no info in the attribute table. Then go the to the other line and catch on the meaning of the line and perform erase between the two named shapefiles. Find them and perform it.How do i do that without making a lot of conditions and line by line: like

directory=input('Insert dir of the main folder')

    If 'INT' in firstline of txt in directory and ASTOM and EAS:
    perform interesection between those.

i don't think the only way is to create hundreads of if conditions that will probably never be valid. How should i do it?

THE LAST VERSION

import os
import fiona
rootdir = r'C:\Users\user\Desktop\a' # path to your root directory you walk
sfiles = {} # a dictionary with all the .shp files
for entry in os.listdir(rootdir):
    dirpath = os.path.join(rootdir, entry)
    if os.path.isdir(dirpath):
        for file in os.listdir(dirpath): # Get all files in the subdirectories
            if file.endswith('.shp'): # If it's an .shp.
                filepath = os.path.join(dirpath, file)
                sfiles[filepath[:-4]] = fiona.open(filepath)

action_dict = {'INT': 'intersection', 'ER': 'difference'}

directory=input('Insert dir of the main folder')
with open(input()) as txtfile: #insert directory of txt
    x = txtfile.readlines()
for line in x:
    action, shape1, shape2 = line.split('_')  # here line is ER_ASTOM_ASTIK or whatever line in your txt file
    if shape1 in sfiles and shape2 in sfiles:
        gpd.overlay(sfiles[shape1], sfiles[shape2], how=action_dict[action])

it doesn't probably perform the intersection or difference. how can i store its results in a variable because i have to do further process with the outcomes of every operation of every line.


texts: How to determine layers affected by a user-drawn polygon

texts:
Scenario:

I have a set of layers, each covering one country. In my OpenLayers client, I enable users to randomly draw a polygon. A user is then able to perform analysis based on the drawn polygon. Analysis include getting the population covered by the polygon, population density of the area covered by the polygon, etc.

This separation of data into separate layers is for performance reasons, as there are millions of grid record data. Putting all this data in one layer is likely to degrade performance considerably.

The Problem:

This works perfectly within a single country, which is represented by a single layer in GeoServer. However, I want to enable users to draw a polygon that covers sections in more than one country. This means that to perform the analysis, data must be extracted from all the affected layers. I'm looking for a pure OpenLayers+GeoServer approach of achieving this. The features extracted are then passed to a GeoServer WPS process.

My attempt so far:

I've created one layer group that includes all the layers for all the countries. I wish to load as minimal data as possible to the client. Which means I can't use WMSGetFeatureInfo, or the WFS version that retrieves information based on a WMS layer. Which means that the layer group approach may not be the appropriate one.

All my processing is done in WPS processes. It's only the result that is sent to the OpenLayers client.

Any suggestions will be appreciated very much.

My environment: GeoServer 2.1.3, PostGIS 1.5, OpenLayers 2.11.


texts: VSIMalloc3() and UnifiedSrcDensity mask errors when running gdal_warp in OSGEO4w shell

texts:
This is related to gdal.py as part of bat script in windows 7

I am getting

  ERROR 2: VSIMalloc3(): Out of memory allocating 854522136 bytes.

  ERROR 2: Out of memory allocating 119813400 bytes for
  UnifiedSrcDensity mask.

I am running a XEON 3.4GHz. Windows 7, 64-bit with 16GB RAM. 930GB on C and 850gb on D (this is free space off 1TB drives) and network drives with 1.72TB free.

The command I am running in a bat file is

gdalwarp -multi --config GDAL_CACHEMAX 500 -wm 500 -r near -co "COMPRESS=DEFLATE" -co "TILED=YES" -co "BIGTIFF=YES" -srcnodata 255 -dstnodata 255 -overwrite %2\%3\%3_translated.tif %2\%3\%3.tif > %2\%3\2_%3_warp_log.txt

where %1=input tif, %2=output directory, %3=file name

In this case the input file is

  Driver: GTiff/GeoTIFF Files:
  P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\4Band_translated_pyramid\z_zones\z_zones_translated.tif
  Size is 411096, 298318 Coordinate System is: PROJCS["WGS 84 /
  Pseudo-Mercator",
      GEOGCS["WGS 84",
          DATUM["WGS_1984",
              SPHEROID["WGS 84",6378137,298.257223563,
                  AUTHORITY["EPSG","7030"]],
              AUTHORITY["EPSG","6326"]],
          PRIMEM["Greenwich",0],
          UNIT["degree",0.0174532925199433],
          AUTHORITY["EPSG","4326"]],
      PROJECTION["Mercator_1SP"],
      PARAMETER["central_meridian",0],
      PARAMETER["scale_factor",1],
      PARAMETER["false_easting",0],
      PARAMETER["false_northing",0],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]],
      EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
      AUTHORITY["EPSG","3857"]] Origin = (15825753.321902087000000,-1808723.335397920600000) Pixel Size =
  (1.000002861523879,-1.000001022298165) Metadata:   AREA_OR_POINT=Area
  DataType=Thematic Image Structure Metadata:   INTERLEAVE=PIXEL Corner
  Coordinates: Upper Left  (15825753.322,-1808723.335) (142d 9'54.58"E,
  16d 2' 4.35"S) Lower Left  (15825753.322,-2107041.640) (142d
  9'54.58"E, 18d35'33.74"S) Upper Right (16236850.498,-1808723.335)
  (145d51'29.19"E, 16d 2' 4.35"S) Lower Right
  (16236850.498,-2107041.640) (145d51'29.19"E, 18d35'33.74"S) Center
  (16031301.910,-1957882.488) (144d 0'41.89"E, 17d19' 5.07"S) Band 1
  Block=256x256 Type=Byte, ColorInterp=Red   Mask Flags: PER_DATASET
  ALPHA  Band 2 Block=256x256 Type=Byte, ColorInterp=Green   Mask Flags:
  PER_DATASET ALPHA  Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA  Band 4 Block=256x256 Type=Byte,
  ColorInterp=Alpha

The pipe of the gdal_warp command gives

  Creating output file that is 411096P x 298318L. Processing input file
  P:\2012\183_TownPlanning_Symbology\Working\Raster_Layer_Creation\4Band_translated_pyramid\z_zones\z_zones_translated.tif.
  Using band 4 of source image as alpha. Using band 4 of destination
  image as alpha. 0

so the process is failing before even 10% is complete and fails about 1h after processing starts.

I am just setting up a linux box to try this out on but is there something I can tweak in windows to get this to work?

best,

--- UPDATE 1

Part of the issue is with the cache and wm switches. Removing these helps but the larger process still crash. http://trac.osgeo.org/gdal/wiki/UserDocs/GdalWarp


texts: How do I snip/split polylines at chainages?

texts:
What is the best way to snip/split lines at chainages?

E.g I have a csv with the uniq ID of the polyline item and a start and end chainage. I want to extract the lines between the chainages stated.

I have access to MapInfo Professional, QGIS and ArcGIS for Desktop, and I am also OK at using MapBasic.


texts: ArcObjects SDK for .NET Framework with Annual Educational Advanced (ArcInfo) Single Use License

texts:
I usually use my organization's ESRI account when at work, but I am now working at my home Desktop PC.  On my home PC, I have an ArcGIS 10.6 annual educational license that I do not want to sacrifice.  It is an Advanced (ArcInfo) Single Use license.  I also have an ESRI Developer subscription through my organization, so I am able to download the ArcObjects SDK for .NET Framework to my home PC.  I also have Visual Studio 2017 Professional on my home PC, which I downloaded BEFORE downloading the ArcObjects SDK for .NET Framework.  When I downloaded the SDK, I ensured that I downloaded the bit for Visual Studio 2017 Professional.  The immediate problem is that Visual Studio is not aware of my ArcObjects SDK.  When I attempt to add a new project, the ArcGIS menu is not on the left side of the New Project dialog box.

Is anybody able to provide a solution, help me troubleshoot or confirm that the SDK will not work with my educational license or maybe it is my Advanced Single Use license that the SDK does not like?

I have read this:  Is the ArcObjects SDK for .NET available with an ArcGIS for Home Use or free trial license?


texts:
Is there a way to to something like ArcGIS provides with the Network Analyst with open source GIS tools?

What I am exactly looking for is:
I am having a roadnetwork and in this roadnetwork I have several points. I want to check if any of these points is further apart from each other than a certain threshold (lets say 1km).

Is there a way to do that with freely available GIS?

I am aware of this
https://anitagraser.com/2011/05/13/catchment-areas-with-pgrouting-driving_distance/
But that just gives me the distance of one point.


texts: Hyperlinking attribute in ArcGIS standalone tables

texts:
I have ArcGIS Desktop 10.0 and have created a 1:m relationship in a File GDB. The 1:m relationship joins a table containing a URL to a point geometry. Some points have more than one file hence the 1:m join. After successfully creating the relationship, I wanted to use the field in the table to open the Image. It seems however, that ArcGIS stand-alone tables cannot use hyperlinks.
Is there a workaround to this or can ArcGIS not use hyperlinks in stand-alone tables?


texts:
I am using QGIS to create maps using CSV files. I want to create a web application where the base map would be first displayed.

I want the users to have dropdowns lists which would be populated with the table attribute values. When the user selects a value from the dropdown a new query would be sent to server and a new layer would then be created on the base map which would be sent back to the client.

Basically I want to learn how do I:

Fetch table attribute values dynamically to web page from QGIS map/shapefile.
Populate those values in respective dropdown lists.
Fire a query to the server where the base map/shapefile resides which the generates a new layer and sends in back to client.

How do I approach this?

I want to code the web app using JavaScript and server using Python.


texts: gdalmerge produces file that is much too large

texts:
I'm trying to merge 60 geotiff files into one huge file in order to be able to product tiles from that resulting huge geotiff file by using gdal2tiles.

The merge process worked, but for some reason it produces a file that has a size of 78GB instead of the expected 19GB.  Not quite sure what went wrong there.

The reason for trying that is because tiling the 60 geotiffs one by one gives me white spaces on the map border.


texts: Navigating ArcObjects Object Model Diagrams (OMDs)?

texts:
How does everyone get round navigating the different Object Model Diagrams (OMDs)?

I can work with single interfaces no problem, but when i have to use multiple interfaces I struggle to get my head round go through different interfaces, as I am never sure which one to use.

So what is everyone's tips that they could share?


texts: Is it possible to export boundaries from Google Map marker?

texts: I primarily have a  doubt that is it possible to export a marked boundaries from google map marker ?if possible kindly tell how to to do this?


texts: Python increasing the filesize of my .mxd

texts:
I have a geoprocessing script that runs the following..

arcpy.AddMessage("TEST...")

Yes just one line of code.

A blank .mxd starts at 391kb for me as I run my script and save the filesize increases.

391 kb
431 kb
451 kb
471 kb
and so on...

A geoprocessing python script seems to increase the filesize by 20-30 kb each time.

I have run this test on an actual script that does something and increases filesize by 40-60 kb...

I have also run on a blank script that does absolutly nothing other than it runs and still increases filesize by 20-30kb...

From what I could find ESRI Tech

  file size increase may be caused by a large amount of geoprocessing results stored in the Results tab in ArcToolbox in ArcCatalog.

Following the steps in the article I was able to decrease the filesize from 589 kb to 293 kb by saving it to a new .mxd

Your thoughts, experiences with this? Has anyone found a way to auto-delete the results tab from a python script so that I don't have to save a new .mxd each time?

UPDATE:

@AHIGH
I have tried playing around with the Geoprocessing options as suggested by AHIGH. Does not seem to stop the filesize from increasing.

@artwork21
If I go to geoprocessing > results and delete the results this WILL decrease the filesize. Is there anyway to automate this process so that the results are auto-deleted after the script has run?

ESRI: Using geoprocessing options to control tool execution

  Any tool results still in the ArcMap table of contents will be saved to disk, regardless of the temporary results setting. This is because geoprocessing has to assume that by saving a map document, you want everything in the ArcMap table of contents saved

I want everything to save except for the results that bloat my filesize...

  If you don't want to save the map document, but you do want to save the data in an individual layer, right-click the layer and click Data > Make Permanent.

Not really a practical solution...

In 10.1 there are new features that would allow me to fix my problem. Guess will have to try to push my company to upgrade to 10.1 ASAP.

Python Tool

WORKAROUND SOLUTION:

Under Geoprocessing > Geoprocessing Options > Results Management > "Keep results younger than:" > Set to Never Save.

This will stop the filesize from increasing.

Alternatively, you can try  to write your scripts so that the run independently of Arcmap (ie. Without it open).

You can write more efficient code, that runs faster. Like a script that runs multiple mxd's at once.


texts:
I am doing something similar to a previous post on stack exchange regarding exporting data from GEE to csv. Once I have the .geo column, how do I work with it in R? I would like to create a shapefile as the final output.
Considering an exported csv as per the previous link called "test" where I have deleted the system:index column that is automatically generated.
test<-{data.frame(B4=c(0.2043, 0.2043, 0.2080), B5=c(0.2573, 0.2573, 0.2653), .geo=c('{"type":"Point","coordinates":[2.6631005755365273,50.90550594166401]}','{"type":"Point","coordinates":[2.6631904070649393,50.90550594166401]}',
'{"type":"Point","coordinates":[2.663280238593351,50.90550594166401]}'))}

df = st_as_sf(test, wkt = ".geo", crs = 4326)

I get this error
Error in CPL_sfc_from_wkt(x) : OGR error
I also tried this to modify the geometry structure
geonew= st_as_sfc(structure(as.character(test$.geo), class ="WKB"),EWKB=TRUE)

I got this error.
Error in CPL_hex_to_raw(vapply(x, skip0x, USE.NAMES = FALSE, "")) :
char2int: unrecognized character in hex string


texts: Extract raster values from point using GDAL

texts:
I want to extract raster value from point shapefile geometry. Raster size is 45924 x 61671 and shapefile has 11.000.000 points. I'm using extract function from Raster package in R but is not very fast.

I'd like avoid non-terminal methods like ArcGIS or QGIS and prefer some terminal options like gdal_xx.py mode at prompt.

As a result I need an object with 11.000.000 length in any format (CSV, txt, shp, etc).

I need repeat this procedure over several similar rasters.

Here's a nice example but I don't know how to get or manipulate the resulting object and I'm not sure about their efficiency.


texts: File naming for Enterprise GIS database?

texts:
I'm designing the guidelines for an Enterprise GIS database.  Part of this database will include customized geoprocessing tools using a variety of software.

I'm interested in knowing whether anyone can foresee a problem with including '_' as a separator for filenames if GRASS GIS is being used as the main geo-processing software?

Most of the tools will be coded in python, SAGA GIS, and R stats will also likely be used as background processess.

Thanks for any advice!


texts: Make all pixel in raster have the same value

texts:
I have a raster file with two categories: category 12 (total count=350) and category 14 (total county = 500). The attribute table looks like this:

rowid Value Count
 0     12     350
 1     14     500

Is there any way I can give all pixel a value of 1?


texts: Make this the only selectable layer (greyed out)

texts:
Someone else has tried Creating shortcut (button) for "Make this the only selectable layer" in Python 2.6? but, as pointed out there, the answer doesn't work.

When adding "Make this the only selectable layer" to the ArcGIS 10.3 menu, using the Customize tool, it adds the button to the main screen, but it is permanently greyed out.

Doing it the old-fashioned way (Table of Contents -> Layer -> Selection -> Make this the only selectable layer) works just fine, but it is a bit slow.

Has anyone had/resolved this problem?


texts: Using ArcPy module in python shell that is not with ArcGIS?

texts:
I have ArcGIS 10.2 installed on my machine. When ArcGIS was installed, python 2.7.3(32 bits) was also installed. I need python 2.7.5 64bits for some gdal functions. So I installed python 2.7.5 64bits stand alone. Now there are two python installed with two python shells. One (2.7.3 32bits) can be started from Start---->All Programs--->Python27--->IDLE(GUI). The other(2.7.5 64bits) is from Start---->All Programs--->ArcGIS--->Python27--->IDLE(GUI).

Till now I have my work done in python2.7.5 64bits. Now I want to integrate a model I built in Modelbuilder ArcGIS 10.2. I export this model to python script so that I can see the code behind. I need to import ArcPy module to use ArcPy functions. I tried several ways but no good luck.

Is there any way that I can import arcpy module in this python shell (2.7.5 64bits, not the one installed with ArcGIS)?

Should I install arcpy package first, like what I did for other module? But I didn't see where to download it.

Finally it is resolved. I followed Is 64-bit background geoprocessing native to ArcGIS 10.2? and asked our IT staff to download the package and installed for me. Now my python (came with ArcGIS) is 2.7.5 64bit. I also uninstalled the stand alone python. Re-install all the modules.


texts: Outputting layer as table at end of process

texts:
I am in the process of migrating some QGIS 2.x processing toolbox script to QGIS 3.x.

I have a collection of scripts that helps in the seismic 3D planning and geometry processing.

One of the scripts is partially working, that is it generates the output layer as expected. For the version 2.x framework I had defined:

##CSV_for_SPS_Generation=output table
##Generated_Lines_Layer=output vector

That CSV is, in reality, a log of the run. It is not meant to be a layer in itself. When running the script the file was generated and automagically loaded as a table without geometry, even as a temporary file.

I was not able to find a solution that could define a file, temporary or otherwise, that would be loaded at the end of the process as a table.


texts: 'Error: table was not found' when running Locate Regions tool using modified raster

texts:
Using ArcGIS 10.7 desktop
I want to run the Locate Regions tool on a land use raster (32 bit signed, with values 0 - 95) in order to identify cell groupings representing a certain land area size. When I run this raster the tool returns an output. But I only want to perform the analysis on upper cell values between 70 - 95. So I used the Extract by attribute tool to get a new raster only containing these values.
I get the following error when I input the 70 - 95 value raster:
ERROR 999999: Error executing function.
The table was not found. [VAT_habitat_suitability_70to95]
No spatial reference exists.
I can't figure out why I get this error. The raster and its table appear to be the same as the original, also 32 bit signed integer and field types are same but just truncated records with the corresponding raster cells. Obviously there is a table and this raster has the same properties, coordinate system, etc.


texts: Reprojecting Points, Shapes and Rasters using QGIS?

texts:
I have a project file that had originally started out with no care to the actual coordinates (we were just using an aerial and displaying the shapes and points on said aerial). Now, I am trying to convert my old coworkers issues, and one was actually using a geographic coordinate system. I have control points to be able to properly scale each layer down and place it in the proper location coordinate-wise, but is there any way of doing this with QGIS? Or should I just get all of the maps redrawn with the correct coordinates in a new file?

An example would be point A - which is currently as 345, -1329 when it should be at -104.671892198, 38.3189727989 (and the transformation should be the exact same for other points, shapes and raster file).


texts: How do you get the nodes of an area in Overpass API in the "right" order?

texts:
I use the Overpass API to get the boundary of an OSM relation. I then use the result to display the boundary on a leaflet map. My problem is however, that the nodes are not in the "right order".

Leaflet of course expects to get a polygon in which the nodes are in the order in which they should be connected afterwards, like in the following picture:

The Overpass API however just gives me the nodes in any order, like for example this:

This problem results in the following output:

Is there some way I can get Overpass to give me the nodes in the "right order"?
I guess that this is the only alternative, as I guess leaflet isn't able to reorder the nodes on its own (it can't really know which order is right).

It has to be possible to get the data in an automatic way. I only need to generate this data once (or maybe once a year) but for over 100 different areas.

Answer to scai's comment:

After closer inspection of the data, not taking the nodes in the order they appear in the result but taking them in the order they appear in the ways seams like the right approach. However, it is still a little more complicated.

There are two problems: First, who guarantees that the ways belonging to a relation are in the right order (all the boundaries I just queried seemed to be ordered the right way, but I have already found a lot of relations where the ways were not ordered in the right way).

Example: Let's assume I get the following data, and the relation consists of

{
    red,
    green,
    blue,
    yellow
}

Secondly, the ways don't all point in the same direction. This problem does definitely happen.

In the example below, the blue way points in the opposite direction than the other ways. Taking the nodes in the order they appear in the ways won't give me the right result.

A possible solution to all this would be the following algorithm: Start with the first way, add all nodes to an array, then search for a new way whose first or last node has the same id as the last added node, repeat until arriving at the first node ever added. But isn't there a better way to do this?


texts: What is a "CFG" application in GIS terminology probably?

texts:
I come across a 2008 tech presentation PDF on "GIS-Enabled Modeling and Simulation (GEMS) presented by ESRI (i think for the US Army).
On one of the page titled "Challenges to be overcomed", there is this point:

> 3D Visualisation:
>    GIS data is inherently 2D
>      Fine for CGF Applications
>    ArcGlobe
>    Procedural terrain visualization approaches.

I believe this point is saying that 2D visualisation is OK for CFG applications, but for their puporse 3D should be used.

Anyone can tell what "CGF Applications" mean?

[edit: deleted second question]

Thanks for any help.

PS: I know it's ESRI sales team job to answer these questions, but I am just a student doing research on terrain rendering. Don't think they will entertain me.

EDIT: typo corrected from CFG to CGF.


texts: Animate and export video over a time series in Google Earth Engine

texts:
I am trying to map a time series of Albedo, and then export it to a video. The mapping works fine, but I am getting the following error on the export:

  Error: ImageCollection must have 3 or 4 bands.

But there is only one band in each of the BlackSkyAlbedo features (images).

The code is included below:

//Map over the dataset
var dataset = ee.ImageCollection('MODIS/006/MCD43A3')
              .filter(ee.Filter.date('2017-10-01', '2018-05-01'));
var blackSkyAlbedo = dataset.select('Albedo_BSA_Band1');
var blackSkyAlbedoVis = {
  min: 0.0,
  max: 400.0,
};
Map.setCenter(-111.475, 45.23, 6);
Map.addLayer(blackSkyAlbedo, blackSkyAlbedoVis, 'Black-Sky Albedo');

print(blackSkyAlbedo);
var polygon = ee.Geometry.Rectangle([-118.79,42.50, -106.0,47]);

// Export
Export.video.toDrive({
  collection: blackSkyAlbedo,
  description: 'blackSkyAlbedo',
  dimensions: 720,
  framesPerSecond: 12,
  region: polygon,
});


texts: Altering existing DEM to approximate its ancient equivalent?

texts:
I'd like to modify a DEM to approximate what an area looked like 3000 years ago, which primarily involves adjusting coastlines based on scattered core sample data.  It doesn't have to be scientifically accurate, although that would be ideal, I'm mostly hoping to achieve this for the aesthetics of an eventual visualization.

I don't really want to go this route, but all I can come up with is photoshop.


texts:
When querying OSM through OverpassTurbo, I notice some locations are ploted as polygons (closed ways) instead of nodes, which makes ill-suited for displaying them in eg. Google My Maps:

WAYS
<Placemark><Polygon><outerBoundaryIs><LinearRing><coordinates>1.2576092,49.2925736 … 1.2576092,49.2925736</coordinates></LinearRing></outerBoundaryIs></Polygon></Placemark>

NODE
<Placemark><Point><coordinates>3.6255999,48.3479532</coordinates></Point></Placemark>

Is there a simple way to turn all polygons into nodes, either with a GUI application like JOSM or a Python script?


texts: How can I call gdal_translate from Java code?

texts:
How can I execute gdal_translate in Java? For example, will something like this work?

gdal_translate "PG:host=localhost dbname='xyzy' schema=public' table='xyz'" "d:\xyx\ab.tiff"


texts: Can't create an ESRI File Geodatabase on one specific drive

texts:
I can create a ESRI file geodatabase in ArcCatalog 10.2 and 10.4 in C:/ drive and in 2 network drives but on the 3rd network drive when you try to create a file geodatabase (fgdb) by either right clicking and selecting new>file GDB, or python script or by using FME. I get an error saying I don't have permission to create a FGDB and it disconnects the drive, not only in ArcCatalog but also if I try to access the drive through Windows Explorer. The only way to re-map the drive is to log out and log back in again.

Upon reopening you can see a new fgdb but it won't open stating it is corrupt or created with a different version to the current version.

It gets stranger because I can create a fgdb using FME's File Geodb-Open API.
I have checked with the network administrators and I have full read/write/execute on the drive and I have no problems creating a personal database .mdb file, shape files or executing python scripts. The administrators have said the drive is fine and there should be no permissions issues and ArcCatalog and ESRI desktop/arcpro all work fine so it doesn't appear to be a software issue either.

There is a slight difference in the size of the corrupted fgdb 115kb compared to 122kb for a normal empty and this appears to be down to the a00000004gdbtable and a00000004.gdbtablx.

How can I fix this problem?

In the short term I can create a fgdb in C drive and copy it across but that's not a preferred long term option. And I would prefer not to have to rebuild my machine.

This question has also been posted on ESRI's geonet website.