texts: SAGA Raster Calculator processing script give error "name 'Rasterout' is not defined"

texts:
I am trying to build a custom script using things I have learned from this post: Number inputs as raster calculator variables in QGIS Modeler

This is my script so far:

##Maechtigkeit unterhalb einer Hoehe=name
##Input_Raster=raster
##Hoehe = string 28.65
##Rasterout = output raster

outputs_SAGARASTERCALCULATOR_1=processing.runalg('saga:rastercalculator', Input_Raster,[],'a - Hoehe',0,False,7,Rasterout)

But I always get the error

  name 'Rasterout' is not defined See log for more details

When omitting the Rasterout parameter and type None instead, it seems to run through but doesn't give me the temporary layer I need? So nothing happens really.

Is there another option for the last parameter? I tried consulting the SAGA documentation but it doesn't yield the information I need.

Also I need that output raster to do more raster calculations on.

How would I be writing that? Can I just rewrite the line:

outputs_SAGARASTERCALCULATOR_1=processing.runalg('saga:rastercalculator', Input_Raster,[],'a - Hoehe',0,False,7,Rasterout)

and use Rasterout instead of Input_Raster? And define a new Output variable'?


texts:
I have two polygons, I want to get the area which does not intersect. I looked at st_disjoint, it returns a boolean, I need a geometry. How do I accomplish this?

Below is my code which gets the intersection area:

INSERT INTO result (id,geom)
(
select id, geom
from way group by id having (
    st_intersection(ST_expand((SELECT ST_Collect(the_geom) FROM vertices
    WHERE id IN( 494733,220418,183786,234791,494733)), +0.1),
    ST_expand((SELECT ST_Collect(the_geom) FROM vertices
    WHERE id IN( 494733,220418,183786,234791,494733)), -0.1))
)
);

My Tables structure

Way Table

Vertices Table


texts: Collecting and analyzing personal movement logs?

texts:
I'm looking for a program that logs what specific locations I've been to & how long I've been there, etc. Is this possible on QGIS?

For example, if I took my laptop to the store today at 1pm (A1) for 30 mins & then went to the gym (B1) for 40 mins at 3pm; the log should look like this:

A1:
Arrived at 1pm and stayed until 1:30pm on 16/03/2019? Yes.

B1:
Arrived at 3pm and stayed until 3:40pm on 16/03/2019? Yes.

This is obviously a very simplified version of what I want. I just want some simple data to work with for the time being.


texts:
I have trouble scheduling a backup in mysql periodically, most of the tools do not recognize the geometric field, for example navicat which says

  Can not get geometry object from data you send to the GEOMETRY field

others simply do not create the table or do not fill the field with the information, someone knows a solution that is free and effective for this problem


texts: Populate markers in leaflet with tweets and tweet location using python

texts:
I am iterating over twitter using python, the user inputs a topic to search for on twitter and then the location to view on the map in leaflet. The leaflet map is created automatically in the browser. How do I get the marker to show up at the tweet location assuming there is coordinates associated with the tweet which in a lot of cases does not happen. Below is the code I think you need to look into this. Suggestions? Using folium, webbrowser, os, geopy, TwitterSearch modules here.

for tweet in ts.search_tweets_iterable(tso):
        tweetText = ( '@%s tweeted: %s' % ( tweet['user']['screen_name'], tweet['text'] ) )
        print tweetText
        if tweet['place'] is not None:
            (lat, lng) = geo(tweet['place']['full_name'])
            print 'And they said it from (' + str(lat) +', ' +str(lng)+')'
        else:
            print "And their place wasn't specified..."
Here is the piece for the marker using the folium module

for tweetLocation in ts.search_tweets_iterable(tso):
    folium.Marker([lat, lng], popup='tweetText').add_to(mapLeafletPython)

I originally had the for tweetLocation loop in the for tweet in ts.search...(tso): part but it wouldn't work. Have them pulled apart doesn't help either.

Here is some documentation for folium:

http://folium.readthedocs.org/en/latest/quickstart.html

also (looking through this now...)

https://pypi.python.org/pypi/folium

Let me know if you have questions, I'd be happy to provide additional details.


texts: When Facebook returns place long/lat, what projection do they use?

texts:
When Facebook Places Search API returns place long/lat, what coordinate system and units of measurement do they use?
According to this question,
In what coordinate system and unit measures the coordinates returned by the Google Maps API are?
This suggests most online organizations use

Coordinate System: WGS84, Web Mercator (Auxiliary Sphere) - EPSG 3857
Planar Units: meters

However, I would like to confirm this.


texts: Tabulate Area in ArcGIS for Desktop not creating correct output?

texts:
I have used the Tabulate Area zonal spatial analyst tool in ArcGIS many times, and am encountering a problem for the first time. I am using a polygon grid (fishnet) as my zone dataset, and tabulating areas of different raster values that fall within the polygon cells, for four different raster datasets (and ultimately comparing the areas of raster values among these different datasets).

For three of my raster datasets, the pixel resolution is finer than the size of the polygon grid cells. The tool works perfectly on these datasets. However, when I run it on the fourth dataset -- for which pixel resolution is equal to that of the polygon cells and the pixels are perfectly aligned with the polygon cells -- the results are completely wrong. The tool should return a table indicating a single raster value associated with each polygon cell. Instead, many cells have multiple values and for these, the sum of areas for each cell are greater than they should be -- 9000 square meters rather than 900 square meters. Yet, when I calculate geometry for the polygon, it returns a value of 900sqm for each cell. I also tried using the zonal statistics tool, and have the same problem.

I have a work-around -- converting the raster to points, then using a spatial join to associate the raster values with the grid cells -- but it is a much more time-intensive process, and I am very curious as to why the Tabulate Area tool will not work in this instance.

Can anyone help  troubleshoot this problem?


texts: Adding non-nullable field to populated feature class

texts:
Having a populated feature class in a geodatabase, a non-nullable field should be added (or a currently nullable field changed to non-nullable).

That current field has currently no null values so I thought it's a straight forward job. With SQL I would use something like:

ALTER TABLE FeatureClassX
        ADD Status INTEGER NOT NULL.
    DEFAULT 0

However, I can not change the properties (non-null) with ArcGIS Pro 2.3.x nor with ArcGIS for Desktop 10.x. Also, using the Add Field and Alter Field tools do not allow me to change this property. Most likely due to the lack of being able to set/config a default value in those tools.

The only workaround I can come up at this point is to create a new feature class with same fields and properties and then copy each record one by one from the original feature class to the new one.

It is not difficult to write a script to do the workaround but I would like to avoid it and hope that I just have missed something in ArcGIS Pro 2.3.x and/or ArcGIS for Desktop 10.x. I am looking for a built-in solution for either of the ESRI products (or using Python) since I am currently lucky to have both installed.

Edit 2019-07-17 0821: Alter Field tool added to non-working solutions after I have tried umbe1987's answer


texts: Converting position data created in R to GPX format?

texts:
Does anyone know if there is a way to convert position data (latitude and longitude) created in R to GPX format to use with navigation software?

I have been searching and havent found an answer.

In R you can bring in GPX formatted data using the readGPS()function in maptools.

I have also found the GPS Visualizer website which will allow me to convert data to GPX.


texts: Selecting multiple polygons - which polygon is top of the merge list?

texts:
If you select multiple polygons, and then try to merge them, ArcGIS creates the dialog box showing the selected polygons (picture below). You can then double click the one you want the others to merge into.

How does ArcGIS decide which polygon comes out top of the list?

I can't work out what criteria are used - the largest polygon, the furthest north, the most recently selected, etc?

If I knew, I could double click on the top one, rather than having to establish which one is which by clicking once to highlight the polygon.


texts:
I have a drone project with 50,000 pictures in 50 layers.
Each layer has 2 to 5 zones.

When I put them all together, we have a lot of overlap and in some places, we have up to 4 pictures of the same area. See the basic example below:

I want to make a mosaic using as little pictures as possible.

How can I use the standard layer management features to manage my imported points and then delete the associated source file or export only the files I need to a new file?

Regards,

Alan C


texts: Deleting first and last vertices using ArcPy of ArcGIS Pro?

texts:
I am trying to delete first and the last vertex of each line.

I created code to delete first and last vertex but it shows me error on line 18:

    Traceback (most recent call last):   File "D:\UKPN\ukpn_CAD\tet\Delete_first_last.py", line 18, in
        insert = arcpy.da.InsertCursor(outFC, "SHAPE@") RuntimeError: cannot open
    'D:\UKPN\ukpn_CAD\ukpn_CAD.gdb\right_no_diss_sample_diss_Sp4' Failed
    to execute (Deletefirstlast).

My code:

import arcpy

arcpy.env.overwriteOutput = True

# set input/output parameters
inFC = arcpy.GetParameterAsText(0)
outFC = arcpy.GetParameterAsText(1)

#output path
fcPath = outFC.rpartition("\\")[0]
fcName = outFC.rpartition("\\")[2]

if arcpy.Exists(outFC):
    arcpy.Delete_management(outFC)

search = arcpy.da.SearchCursor(inFC, "SHAPE@")
insert = arcpy.da.InsertCursor(outFC, "SHAPE@")

for row in search:
    points = [arcpy.Point(point.X, point.Y) for shape in row[0] for point in shape]

    del points[0]
    del points[-1]
    points.append(points[0])

    array = arcpy.Array(points)
    polyline = arcpy.Polyline(array)

    insert.insertRow([polyline])

I am using ArcGIS Pro Basic license 1.3. Can anyone tell me what am I doing wrong?


texts: GDAL changes pixel resolution when transforming from HDF to netCDF

texts:
Using the MODIS product MCD18A1.061 (DSR 1km), tile h15v05-2001001, my goal is to extract the scientific datasets that are contained in an HDF to several netCDFs.
For this, I used the following gdal commands:
gdal_translate -of netCDF -co 'FORMAT=NC4' HDF4_EOS:EOS_GRID:'{input}':MODISRAD:{sds} {output}

gdalwarp -of netCDF -co 'FORMAT=NC4' {input} {output} -s_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -t_srs '{shapeproj}' -srcnodata -1 -dstnodata -1 -overwrite

""gdalwarp -of netCDF -co "FORMAT=NC4" -cutline {shapefile} -csql "SELECT * FROM {shapefilename} WHERE {prop1}='{prop2}'" {input} {output} -srcnodata -1 -dstnodata -1 -overwrite""

According to the following images:

My problem seems to be in the first command (gdal_translate), since when converting from HDF to netCDF the resolution of the pixels decreases. I would like to know why the resolution is changed and how can I prevent this so that the netCDF file has the same resolution as the original (1km).
I tried to use the -tr option but got the same results, only the file size changed.


texts: Can you define a topology rule between two feature classes from different feature datasets?

texts: I work in a project that I want to define topology rules between feature classes from different feature datasets, is that possible ? if no is there a solution that can do that ?


texts: arcpy.intersect_analysis: output shapefiles not saved?

texts:
I am trying to write my first python code to use in ArcGIS 10.1.

I have feature classes in a file geodatabase .gdb, and AOI.shp (area of interest). I list my geodatabase feature classes, and want to intersect my list with AOI.shp.
I want to save my final intersected files.

Here is my code:

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

# Allow to overwrite existing files
arcpy.env.overwriteOutput = True

# Set environment settings
env.workspace = r"C:\Users\...\extent_AOI\disturbaces_for_aoi.gdb"

# identify AOI
aoi = r"C:\Users\...\aoi.gdb\aoi"

# Set output database
out_wksp = r"C:\Users\...\disturbances_by_aoi.gdb"

# list features in the working folder database
fcl = arcpy.ListFeatureClasses()

# print files to check them
print fcl

for fc in fcl:
    arcpy.Intersect_analysis([fc, aoi], out_wksp + fc, "ALL", "", "INPUT")
    print "intersection ok"

My code seems to run through PythonWin.

However, I can't find my output data in my output folder out_wksp ?

How can I save my output data correctly?


texts: ogr2ogr - accomodate spaces in -select properties

texts:
I have a geojson feature collection. Some of its features have spaces in their property names. I am attempting to use the -select argument to limit which properties are included in the output.

Is there a way to pass in property names that have spaces in them? On Ubuntu 18.04.3 I have tried:

ogr2ogr -select '"My Property"' -f 'Esri Shapefile' my-shp-dir my-fc.json

Which results in: Error 1: Field "My not found in source layer

I have also tried not quoting the spaced properties, and putting square brackets around them.

Is there a way to handle spaces in properties?


texts:
Why if I applied a code over a MODIS collection and then I download an image, the image is in WGS projection? The code is the next:

// Simple regression of year versus NDVI.

// Define the start date and position to get images covering TDPS,
// Bolivia, from 2000-2010.
var start = '2000-01-01';
var end = '2008-12-31';

var region = table3

// Filter to Collection
// time of year to avoid seasonal affects, and for each image create the bands
// we will regress on:
// 1. A 1, so the resulting array has a column of ones to capture the offset.
// 2. Fractional year past 2000-01-01.
// 3. SAVI.
// 4. Mask Clouds
var col = ee.ImageCollection("MODIS/MOD13Q1")
  .filterDate(start, end)
  .filter(ee.Filter.calendarRange(1,1,'month'))
  .filterBounds(region)
  .map(function(image) {
    var date = ee.Date(image.get('system:time_start'));
    var yearOffset = date.difference(ee.Date(start), 'year');
    var savi = image.expression(
      '(1 + L) * float(nir - red)/ (nir + red + L)',
      {
        'nir': image.select('sur_refl_b02'),
        'red': image.select('sur_refl_b01'),
        'L': 0.5
      });
    return ee.Image(1).addBands(yearOffset).addBands(savi).toDouble();
  });

print(col)
// Convert to an array. Give the axes names for more readable code.
var array = col.toArray();
var imageAxis = 0;
var bandAxis = 1;

// Slice off the year and savi, and solve for the coefficients.
var x = array.arraySlice(bandAxis, 0, 2);
var y = array.arraySlice(bandAxis, 2);
var fit = x.matrixSolve(y).clip(region);

// Get the coefficient for the year, effectively the slope of the long-term
// SAVI trend.
var slope = fit.arrayGet([1, 0]).clip(region);

Map.addLayer(slope);

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: slope,
  description: 'SAVI_slope_jan_MODIS_00_08',
  scale: 250,
  maxPixels: 1e9,
  region: geometry
});

But if I download an image from that collection it is in MODIS SINUSOIDAL. The code is the next:

//Choose country using GEE Feature Collection
var region = table;

var collection = ee.ImageCollection('MODIS/MOD13Q1')
.filterDate('2000-01-01', '2014-12-31')
.filterBounds(region)
.filter(ee.Filter.calendarRange(1,1,'month'));

print(collection)

// Collect data and filter by total dates
var modis_savi = ee.Image('MODIS/MOD13Q1/MOD13Q1_005_2001_01_01');

print(modis_savi)

var savi = modis_savi.expression(
      '(1 + L) * float(nir - red)/ (nir + red + L)',
      {
        'nir': modis_savi.select('sur_refl_b02'),
        'red': modis_savi.select('sur_refl_b01'),
        'L': 0.5
      });

Map.addLayer(savi.clip(region))

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: savi,
  description: 'MODIS_JAN_1_2001',
  scale: 250,
  maxPixels: 1e9,
  region: table
});

Does anyone know why?


texts: API to reverse geocode latitude, longitude to census tract?

texts:
I have a bunch of latitude,longitude coordinates (mostly city/town centroids) in the United States for which I'd like to find the containing census tract.

For instance,

input => nn.nnn, -nn.nnn
output => 010300

Is there a free or low-cost API out there that can help? (Essentially Yahoo Placefinder, with the additional response of census tract.)

(I understand that many cities include multiple census tracts and I'd expect that some small towns contain only parts of one or more census tracts. I'm comfortable ignoring that issue for now.)


texts:
I'm currently working with OS Code-point polygons and google maps coordinates. I have google maps point coordinates for a list of houses which are located at the center of each full postcode (according to google maps) they correspond to. The code-point polygon OS data has polygons for each full postcode marking their boundaries, but also has a fair few (5%-ish) small polygons, called "vertical streets" - small areas (usually tall buildings) within other post code polygons which have more than one postcode assigned to them. These each have a reference number and separate look-up table which gives the (usually) 2-3 postcodes that are assigned to that building.

How can I try and assign a postcode in buildings that fall in the vertical streets, given that I don't have a full address (only street name) for these properties?

I have checked online and can't find any indication of how this can be dealt with. Google maps postcode data is not the most accurate either - I'm not sure how to deal with this.


texts: How to get elevation profile from .tif dsm on the web map?

texts:
I have my dsm .tiff file, pixels are floating point 32bit values giving the terrain elevation in metres. I need to create a web map where user can draw a path and get elevation profile for that path. Raster has attribute 'GRAY_INDEX' where values are stored.

What are the options to achieve this and is it even possible?

I thought about Geoserver GetFeatureInfo, but I'm not sure how that would work as I can get values for only one point at a time. I thought when user draws a path, I make a number of GetFeatureInfo requests exactly as points in that path, but I assume that would lead to slow performance. Is there a better option?


texts:
My company is interested in using GIS to create sales territories in the pharmaceutical industry. We currently use "alignment" software, but now we're looking for more commonly used and versatile GIS options.

Primarily, we would like the ability to combine zip codes, cities, or counties into territories, and easily drag and drop them into other territories.


texts: Converting Excel Spreadsheet into Polygon Shapefile

texts:
I am trying to convert the Wards in Leeds data in data.gov.uk into a polygon shapefile, so I can clearly see all of the wards and perform calculations.
I have tried using the add spreadsheet layer option, but all I get is a dot on the screen. I have tried the add vector layer option too, but I either end up with another point or nothing at all. I am making sure that the co-ordinate system I'm using is correct.
If you could provide a step-by-step guide of how to convert the excel spreadsheet into a polygon shapefile so I can see all of the wards, I'd be very grateful.
I am working with QGIS 3.14, and you can find the data at: https://data.gov.uk/dataset/ec9e6e2f-d4c0-4d70-9106-c9a4eff49ff3/wards-in-leeds


texts:
I am attempting to use eCognition Developer to segment and classify 5-band UAS imagery. When I go to create a new project in eCog and input my imagery (individually selecting the five .tif files), it is creating duplicate layers. I should only have 5 layers in my project, but it is returning 11 layers. Red, Green, Blue, NIR, and Red Edge layers are all duplicated, and an extra (undefined) band is input with them.

Below is a screenshot of my Image Object Information screen (after a segmentation run) so you can see what layers I am referring to. I have my layers named according to what they represent, and the duplicates are just "Layer 1,2,etc."

Is there any way I can prevent this from happening, so I only have my 5 layers of interest within my project?


texts: OpenLayers Lat Lon overlay animation doesn't work

texts:
I'm new to OpenLayers and trying to experiment with this code snippet. The idea is to show the Los Angeles area and then overlay it with several points defined at certain latitude and longitude. Highlight and animate them to show how these geo-codes are important to the business. I've changed this code to:

Zoom into the Los Angeles area
Instead of the random points, I've defined two arrays of latitudes and longitudes. Then, I want a lat and long to be chosen from these arrays but only correspondingly. So not (lon[0], lat[2]) but only at the same index like (lon[2], lat[2]). And then that point defined by longitude and latitude should be animated.

For these two changes, I altered the original js code at these two places:

view: new ol.View({
          center: [-13180189,4018753],
          zoom: 8

function addRandomFeature() {
        var lat = [34.152568,34.069787,34.21377,33.81761,33.82508,33.824989,33.905536,33.895896,33.921715,34.19184,34.190977,34.021206,34.050864,34.060165,34.07235,34.076094]
        var lon = [118.347831,118.384738,118.227904,118.070374,118.14034,118.144561,118.05531,118.046738,118.098677,118.33581,118.338151,117.976882,117.989626,117.991428,118.381041]

        for (i=0; i< lat.length; i++){
            x = lon[i];
            y = lat[i];
            var geom = new ol.geom.Point(ol.proj.fromLonLat([x, y]));
            var feature = new ol.Feature(geom);
            source.addFeature(feature);
        }
      }

Zooming into Los Angeles works but the second change doesn't show the points defined by the lat long highlighted as the blue spot that gets animated, as was the case in the original code snippet -- albeit defined by random points in the space, not latitude and longitude. A single point somewhere in China is being shown with a red halo. No lat long in the Los Angeles area is being highlighted with the blue circle.


texts: Convert latitude and longitude in shapefile to different scale

texts:
I'm working with a shapefile where the latitude and longitude is large.

Here is it bounding box:

> utah_ed@bbox
                min       max
coords.x1  239109.1  671324.5
coords.x2 4101331.7 4651097.8

However, when getting a Utah map from ggmap, the longitude is between (-113, -110) and the latitude is between (39.5, 42.0).

How do I convert them into the same scale so I can plot over them?


texts: R will not read shapefile / file.exists returns true

texts:
I cannot successfully read a shapefile in R.

I'm using Rstudio version 0.99.902

I have read several articles on readOGR errors that concern reading a shapefile.

My shapefile is error free, as it opens in qGIS. when I check prescence of file, file.exists returns TRUE yet I get a "Cannot open layer" error:

    Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = >>use_iconv, : Cannot open layer

Any advice?


texts:
In a for loop, the code tried to select one dam feature out from a set of dams and create a new shapefile so to perform cost path analysis with it. However, it seems field value of the new shapefile cannot be extracted or used to perform a cost path analysis. The shapefile was selected and created correctly and I was able to see fields when adding data to ArcMap.

Here's the code:

rows = arcpy.SearchCursor(inPtFeatures)
for row in rows:
  gId = row.getValue(keyAttr)

print ("Processing "+str(gId)+"...")
elev = row.getValue(keyAttr2)

  if os.path.exists(procDir+"/"+str(gId)):
    shutil.rmtree(procDir+"/"+str(gId))

os.mkdir(procDir+"/"+str(gId))
currentDir = procDir+"/"+str(gId)
currentPtFeat = currentDir+"/pt_"+str(gId)+".shp"

whereSQL = '"GRAND_ID" = '+str(gId)
arcpy.Select_analysis(inPtFeatures,currentPtFeat,whereSQL)
print ("test2")

method = "EACH_CELL"
destField = "FID"
outCostPath = arcpy.sa.CostPath(currentPtFeat,inDem,inFdir,method,destField)
outCostPath.save(currentRastPath)

and error message:

\Functions.py", line 688,
in CostPath     destination_field)   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py", line 47,
in swapper     result = wrapper(*args, **kwargs)   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 681,
in wrapper     destination_field)   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 498,
in <lambda>     return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 000860: Input raster or feature destination data: is not the type of Composite Geodataset, or does not exist. ERROR 001000: Destination field: Field FID does not exist


texts: How can I transform polylines into points every n metres in PostGIS?

texts:
Question in pretty much the same as this, but I need to use PostGIS.

I have a number of polylines, and want to trasform them into points every n metres along these lines.

edit:
Thanks a lot to Stefan. My final query was a bit different, probably due to older version of PostGIS (1.5), but in the end I got what was needed. Points were combined into a heatmap of roads density.


texts: I am trying to run the test files for flow mapper plugin in QGIS but I run into an error message.


texts: Automatically calculate field with every new entry?

texts:
I have a point feature class that has a "Direction" field (north, northwest, west, etc; linked to a domain with 16 options). I then have a numeric field "Rotation" that controls the rotation of the feature, based on the string in the "Direction" field.

The value of the numeric field "Rotation" is based on a simple python script of "if" conditions to determine the numeric value of "Rotation" based on the strings from "Direction".

For example, if the point has "Direction" field 'North', the rotation field is 0. If it is 'South', it is 180, and so on.

Is there a way to set up the "Rotation" field to update automatically for every new entry in the attribute table?


texts:
I have a .csv file of coordinates such as the following:

-90.19, 38.62, St Louis
-122.42, 37.82, San Francisco
-74.04, 40.68, New York

I want to create a shapefile using DotSpatial that has these points, so that I can overlay them on a shapefile of the United States.  Which functions in DotSpatial would I provide these points to, to generate a shapefile?

I don't use this library, however you can write a similar code:

            string[] csv = new string[] { "-90.19, 38.62, St Louis", "-122.42, 37.82, San Francisco", "-74.04, 40.68, New York" };

            using (FeatureSet features = new FeatureSet(FeatureType.Point))
            {
                features.CoordinateType = CoordinateType.Regular;
                features.Filename = @"c:\temp\test1.shp";
                features.Projection = ProjectionInfo.FromProj4String(KnownCoordinateSystems.Geographic.World.WGS1984.ToProj4String());
                DataColumn dt = features.DataTable.Columns.Add("City", typeof(string));
                int idxColumn = dt.Ordinal;
                foreach (var row in csv)
                {
                    string[] s = row.Split(',');
                    IBasicGeometry g = new Point(Convert.ToDouble(s[0], CultureInfo.InvariantCulture), Convert.ToDouble(s[1], CultureInfo.InvariantCulture));
                    IFeature f = features.AddFeature(g);
                    DataRow dr = f.DataRow;
                    dr.SetField<string>(idxColumn, s[2].Trim());
                }

                features.Save();
            }


texts: CartoDB stripping the address from KML files

texts:
I am trying to import a kml file into CartoDB but everytime I do it seems to strip the address information.

Here is a sample of the KML file:

<Placemark>
            <styleUrl>#0</styleUrl>
            <name>AOJI Educ</name>              <ExtendedData>
                                            <Data name='Country'>
                            <value>Australia</value>
                        </Data>
                                            <Data name='Phone '>
                            <value>02 9264 7171</value>
                        </Data>
                                            <Data name='eMail'>
                            <value>aileen@globee.com</value>
                        </Data>
                                            <Data name='Website'>
                            <value>www.globee.com</value>
                        </Data>
                                            <Data name='Main Contact'>
                            <value>Ailen Wan</value>
                        </Data>
                                </ExtendedData>
                            <address>Suite 8, Level 6, 5961 George Street, Sydney New South Wales 2000 Australia</address>
            <Point>
                <coordinates>151.2058109,-33.8774203,0</coordinates>
            </Point>
        </Placemark>


texts: Appropriate web-server framework for serving slices of an interactive Leaflet vector layer from a large PostGIS table?

texts:
I'm trying to develop a web map with a clickable chloropleth layer of streets in a large city. There are approximately 40k polylines in this layer, producing a 21MB GeoJSON. This is too much data to send to the client at once and an unnecessary amount to map at the city scale. Our data is in PostGIS and we'd like to Leaflet for the front-end. I have investigated the more common Python web-servers but I haven't easily found documentation/tutorials/Q&A about solving something similar with geoDjango/geoAlchemy2.

The answer to Adding/Removing Leaflet GeoJSON layers shows how to load different layers at different zoom levels (and caching them), but doesn't cover a bounding-box filter on data coming from the server. This answer has the rather unhelpful solution of "just load all the data client-side first". The answer to How to add a bounding box filter to this leaflet WFS request? shows how to load data within the bounding box on every map move, but, unlike the first answer, appears to discard and re-draw the layer every time the map is moved.

geojson-vt by Mourner solves the client-side rendering problem (106MB of zipcodes, wow!) but doesn't seem to solve the data-transfer issue. The post does hint that this has traditionally been done server-side.

  The best way to optimize the data for all zoom levels and screens is to cut it into vector tiles. Traditionally, this is done on the server, using tools like Mapnik and PostGIS.

What python web-server framework can load and display slices of a large number of polylines within different zoom levels & bounding boxes?


texts: WKT to some object that will give me longitude and latitude properties

texts:
How, in JavaScript and OpenLayers, can I pass my WKT point string to a method and get back an object that will parse the longitude and latitude into properties for me?

'POINT (XXX YYY)'
to
SomeObject.longitude
and
SomeOjbect.latitude


texts: Extracting feature centroid in shapefile using FME?

texts:
I have 50 shapefiles in point, line and polygon.

Each layer have different code; for example road have code 1120 and tracks have 1121 and streams have code 2110.

I want to calculate centroid of each feature up to 9 figure each for longitude and latitude prefix by layer code.

I want output as 1120790989898329085438.

In above code "1120" is feature code, longitude is "79.0989898" and latitude is "32.9085438".

Is there any way to calculate centroid in FME and output in new field or defined field?


texts: How to reorder field in a file geodatabase without creating a new feature class?

texts:
I have a Script(toolbox) to reorder the fields into a file geodatabse and the output is a new feature class with fields in the order specified, is there any way to replace the existing feature by the new input created?
I'm using ArcGIS 10.2
here is my script :

import os,string,sys

try:
    import arcgisscripting
    gp = arcgisscripting.create()
except:
    import win32com.client
    gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")

FC_In = sys.argv[1]
NewOrder = sys.argv[2]
FC_Out_WS = sys.argv[3]
FC_Out_Name = sys.argv[4]

strFields = gp.listfields(FC_In)
strNewOrder = NewOrder
Layers = [""]

# Get Field Collection into List
Fields = gp.listfields(FC_In)

NewOrder = string.split(strNewOrder,";")
print Fields
print NewOrder

x=0 #NewOrder enumerator
y=0 #Layers enumerator (first)

# Cycles through the input feature class creating field collections
while x<len(NewOrder):
    Fields.Reset()
    Field = Fields.Next()
    while Field:
        FieldName = str(Field.name)
        if x<len(NewOrder) and NewOrder[x].lower() == FieldName.lower():
            Layers[y] = Layers[y] + FieldName + " " + NewOrder[x] + " VISIBLE"
            x=x+1
        else:
            Layers[y] = Layers[y] + FieldName + " " + FieldName + " HIDDEN"

        Field = Fields.Next()

        if Field:
            Layers[y] = Layers[y] + "; "
            FieldName = str(Field.name)
        else:
            print Layers[y]
            Layers.append("")
            y=y+1
            Layers[y]=""

# Use Field collections to create Layers with proper order until all fields have been used
z=0 #Layers enumerator (second)
LayerList = ""
for Layer in Layers:
    if Layer <> "":
        if LayerList <> "":
            LayerList=LayerList + ";"
        LayerList=LayerList + "Layer" + str(z)
        print LayerList
        gp.MakeFeatureLayer(FC_In, "Layer" + str(z),"#","#",Layer)
        gp.addmessage(gp.getmessages()+ "\n")
        z=z+1

#Collect Info about Input Feature Class
dsc = gp.describe(FC_In)
fc_in_stype = dsc.shapetype
fc_in_sr = dsc.spatialreference

# Create new empty Feature class with desired fields
try:
    gp.CreateFeatureClas(FC_Out_WS,FC_Out_Name,fc_in_stype,LayerList,"#","#",fc_in_sr)
    gp.addmessage(gp.getmessages()+ "\n")
except:
    gp.addmessage(gp.getmessages()+ "\n")

# Append Records to new feature class
print FC_Out_WS + "\\" + FC_Out_Name
print FC_In
gp.append_management(FC_In, FC_Out_WS + "\\" + FC_Out_Name)
gp.addmessage(gp.getmessages()+ "\n")

del gp


texts: Natural Breaks Results Difference in Different GIS Analyst Tools

texts:
I calculated natural break clustering in each GIS software and Python Jenkspy library.

The data is population of each city in Europe and this is the result below.

Jenkspy library and QGIS shows somewhat similar results (but not really) while other trials return totally different results.

I am going to try with CartoDB PostGIS Extension though, I am already doubting the result.

the input data is all same and I made 6 classes for division. But I am not understanding why the results are different.

Does anyone know the reasons?

Surprisingly,
QGIS 2.18 and QGIS 3 gave me different results with same data.


texts:
I have some confusion about RoadLink and RoadArea topological relation in INSPIRE (Infrastructure for Spatial Information in the European Community).

If I got this correctly:

RoadLinks are connected at the intersections.
RoadLink must be inside it's RoadArea.
RoadAreas must not intersect.

So my question is: How can I solve topology on the intersection? If I keep the RoadLinks connected then the RoadAreas will intersect. If I trim the RoadAreas (dashed line between two RoadAreas at the picture below) than RoadLink will be outside it's RoadArea.


texts: OWSLib - ValueError: All strings must be XML compatible: Unicode or ASCII

texts:
I try to use OWSLib to search metadata in GeoNetwork, and I have problems when the string contains non-ASCII characters.

geonetwork = 'http://geonetwork-mshe.univ-fcomte.fr:8080/geonetwork/srv/fre/csw'
from owslib.csw import CatalogueServiceWeb
csw = CatalogueServiceWeb(geonetwork)

1) ASCII string. It is OK.

myquery = PropertyIsEqualTo('csw:AnyText', 'africa')
csw.getrecords2(constraints=[myquery], maxrecords=20)
csw.results

  {'matches': 1, 'nextrecord': 0, 'returned': 1}

2) string with accents

myquery = PropertyIsEqualTo('csw:AnyText', 'été')
csw.getrecords2(constraints=[myquery], maxrecords=20)

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/owslib/csw.py", line 355, in getrecords2
    node2.append(flt.setConstraintList(constraints))
  File "/usr/lib/python2.7/dist-packages/owslib/fes.py", line 151, in setConstraintList
    flt = self.setConstraint(constraints[0])
  File "/usr/lib/python2.7/dist-packages/owslib/fes.py", line 120, in setConstraint
    self._root.append(constraint.toXML())
  File "/usr/lib/python2.7/dist-packages/owslib/fes.py", line 280, in toXML
    etree.SubElement(node0, util.nspath_eval('ogc:Literal', namespaces)).text = self.literal
  File "src/lxml/lxml.etree.pyx", line 1031, in lxml.etree._Element.text.__set__ (src/lxml/lxml.etree.c:55347)
  File "src/lxml/apihelpers.pxi", line 711, in lxml.etree._setNodeText (src/lxml/lxml.etree.c:24667)
  File "src/lxml/apihelpers.pxi", line 699, in lxml.etree._createTextNode (src/lxml/lxml.etree.c:24516)
  File "src/lxml/apihelpers.pxi", line 1439, in lxml.etree._utf8 (src/lxml/lxml.etree.c:32441)
ValueError: All strings must be XML compatible: Unicode or ASCII, no NULL bytes or control characters

I tried to use &eacute; instead of é, but without success...

Any ideas?


texts: having a standard value when filling in attribute table

texts:
I am using QGIS to collect field data and I wonder if it is possible to have a standard value apply to some of my attributes unless I say otherwise.

I fill my attribute table with several drop down menues ("name", "species", "quantity"..) but in most cases, some of these attributes are the same over and over again.
EDIT: I use Value Relations to populate the menues..

So I wonder if there is a way to either

a) always have, lets say the first entry of a specific drop down (attribute) pre-selected upon point creation

or

b) always fill the attribute table with a specific value if I select nothing. e.g. if I select nothing specific for the "quantity" attribute, fill the table with "1"


texts: QGIS: connection to Snowflake cloud warehouse to accessing data tables

texts: Is there a way in QGIS to connect to a Snowflake cloud data warehouse using ODBC, or other? I would like to access this data storage to natively create spatial content by accessing the WKT, Geography, etc. Or by creating an X,Y using the lat/longs stored in the table.


texts:
Each day I need to draw a path on a map and add a text (like 4, 5 or 8 min) indicating how long it takes by car from the starting point to the destination (see the figure below).

I thought it would be helpful to create a Shiny App using Leaflet in R (code is shown below).

I make use of addDrawToolbar from the leaflet.extras package to draw the path as can be seen on the map attached.

But I do not know and could not find how to add a text in the same manner as I draw the path. The solution does not strictly need to be in R.

My aim is to create an App for someone who would like to do these kinds of things and at the same time who does not know how to code.

library(shiny)
library(leaflet)
library(leaflet.extras)

ui = fluidPage(
      tags$style(type = "text/css", "#map {height: calc(100vh - 20px)
      !important;}"),
      leafletOutput("map")
      )

server = function(input,output,session){
             output$map = renderLeaflet(
                 leaflet()%>%

         addTiles(urlTemplate = "http://mt0.google.com/vt/lyrs=m&hl=en&x=
              {x}&y={y}&z={z}&s=Ga")%>%

         addMeasure(
              primaryLengthUnit = "kilometers",
              secondaryAreaUnit = FALSE
         )%>%

         addDrawToolbar(
              targetGroup='draw',

              editOptions = editToolbarOptions(selectedPathOptions =
                    selectedPathOptions()),

              polylineOptions = filterNULL(list(shapeOptions =
                    drawShapeOptions(lineJoin = "round", weight = 8))),

              circleOptions = filterNULL(list(shapeOptions =
                    drawShapeOptions(),
                    repeatMode = F,
                    showRadius = T,
                    metric = T,
                    feet = F,
                    nautic = F))) %>%
        setView(lat = 45, lng = 9, zoom = 3) %>%
        addStyleEditor(position = "bottomleft",
                 openOnLeafletDraw = TRUE)
 )
}

 shinyApp(ui,server)


texts: Changing Data Types within the CMV Map Service

texts:
Is it possible to format fields that appear in tables within the map service?

In the table results after a spatial query, date fields are converted to UTC rather than duplicating  the values from the original data (e.g., original data is "1/2/2017"; the table in the CMV shows "1/1/2017  6:00 p.m.").  Export.js was modified with a Date() function so the exported file has the correct date.  Also, a field that is a multi-digit integer value is appearing as a floating point with decimal places and commas (e.g., value:  123456,  table field: 123,456.000). This value is a code and isn't a decimal value.


texts: Pgr_WithPoints : POIs in different tables

texts:
Is it possible to calculate the minimal distance between 2 points outside of routing network using pgr_withPoints, knowing that these points are from different tables?
I tried to group them in one table and it did work, but for efficiency I need to keep them separated in two tables


texts: Merging two DEMs in different pixel bits?

texts:
I have 2 different Digital Elevation Modeling (DEM). One is the state DEM and has 30m and has 16 pixel bit. I have a couple of smaller DEM that has 30m but they have 32 pixel bit.

I know how to merge them, but in my mind, do the Pixel bits have to be in the same pixel bits for example, 30m 16 pixel bits equal to 30m 16 pixel bits, correct ?

Reason for asking this question was because when I am working on a watershed and clipped state 30m 16 bit in half and the other smaller DEM has 30m 32 bits and merge them.

I tried to find this on the internet and found nothing.


texts: Programatically retrieve units from a map's SpatialReference

texts:
I'm trying to retrieve the units (like feet or meters) of the current spatial reference from an ArcGIS map, using the JavaScript API.

I have the map's spatial reference object, with it's well-known ID.  I'm guessing there would be a lookup function or something I could use to get the units for that spatial reference, but I can't find one.

Does such a function exist?


texts: Viewing list of data sets that have been edited under specific ArcSDE Version?

texts:
I'm wondering if it's possible to list all of the data sets (feature classes, tables, etc.) that have been edited under a particular version in ArcSDE.

A tool that does this would be ideal (I see there used to be one by NASA, but it doesn't seem to be working for me in 10.3), but I am also open to any other suggestions, including accessing SDE directly.

And just to be clear, I am not looking to compare changes made between a version and the default. I already know how to do that using the "Version Changes" button on the Versioning toolbar in ArcMap.


texts:
I'm having issues with ArcGIS Pro (AGP) crashing when I run this ArcPy code directly into ArcPro Python window. I previously wrote the three lines separately but I decided to try to copy python code from Geoproccessing history. However the outcome is the same where Pro Crashes. Specifically, the trailing dots freeze, and the program freezes and the below image pops-up. There is no error code thrown, just a Windows Program 'has stopped working' as seen below

Does it have something to do with how I'm storing or not storing the data?

The code is supposed to select some attributes, copy it to a new layer, and dissolve that new layer.  They seem to work individually when run.

import arcpy

arcpy.env.workspace = 'c:\data\python.gdb'

arcpy.management.SelectLayerByAttribute("SHAPE_START", "NEW_SELECTION", "LEASE_NAME LIKE 'ABC%'", None)

arcpy.management.CopyFeatures("SHAPE_START", r"C:\Data\GIS\Project\Python\Python.gdb\INPUT_SHAPE_SELECTION", None, None, None, None)

arcpy.management.Dissolve("INPUT_SHAPE_SELECTION", r"C:\Data\GIS\Project\Python\Python.gdb\FINAL_OUTPUT", None, None, "MULTI_PART", "DISSOLVE_LINES")

I am using ArcGIS Pro (AGP) version 2.2.4 on Windows 10 using the ArcGIS Pro Python Window


texts: Python error can't find PyQt.Qtcore, when launching QGIS

texts:
I installed QGIS in Linux and it was working fine. I try to launch the program again and I get the following Python error:

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

Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 219, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/natalia/.qgis2/python/plugins/MetaSearch/__init__.py", line 29, in classFactory
    from MetaSearch.plugin import MetaSearchPlugin
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 478, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/natalia/.qgis2/python/plugins/MetaSearch/plugin.py", line 28, in
    from qgis.PyQt.QtCore import QCoreApplication
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 478, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named PyQt.QtCore

However, I do have installed PyQt.QtCore. I think the problem is with my path and QGIS doesn't know how to look for the Python package.

I am really new to this, I am not sure how to fix it.  The way that my path looks is:

PATH=/home/usr/.qgis2/python:$PATH
PATH=/home/usr/share/qgis/python:$PATH


texts:
While importing NDVI dataset in ArcGIS should we also import the VI_Quality dataset ?

If yes then how should we handle it ?

and, if not, then for what it is there ?


texts:
I am currently stuck on a problem surrounding ArcGIS Explorer and comboboxes

What I am trying to do is use 3 comboboxes to do the following:

View all features contained within ArcExplorer's table of Contents
View the selected layers containing field names
View unique attributes for the fields

All of these would go in 3 seperate comboboxes.

Would anyone have any suggestions on working through this?


texts: GRASS New/Open Mapset not available ( QGIS 2.18.7) on Mac OS

texts:
I have QGIS 2.18 and the option to create a New Mapset is unavailable to me. GRASS is installed but when I go to the Plug-Ins>GRASS only Open GRASS Tools and GRASS Options are available to me. These are also the only icons in the GRASS toolbar.
How can I work with GRASS?

Update, I also notice this ('Currently selected GRASS installation is not valid'):


texts: Convert Shapefile from EPSG:5320 to 4326

texts:
I am trying to need Electoral District boundary data available at the below-mentioned link to SQL Server.

https://www.elections.on.ca/en/voting-in-ontario/electoral-district-shapefiles/open-use-data-product-licence-agreement/download-shapefiles.html

I downloaded the 2018 General Election files and checked the CRS through QGIS 2.10.1 as well as the http://prj2epsg.org/search website. QGIS shows the CRS as EPSG:5320, NAD 83. Whereas the prj2epsg.org shows a number (4540!) of EPSG matches!

Anyways, I tried to change the EPSG:5320 to EPSG:4326 using ogr2ogr as well as QGIS.

For changing the ESPG through QGIS, I followed the steps through 10 as mentioned in the link below.

https://gist.github.com/peter-rose/838a8fe3a725261527334f04cfc3dc00

I also tried to convert the EPSG:5320 to EPSG:4326 using ogr2ogr as well.

ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:4326 -s_srs EPSG:5320 ELECTORAL_DISTRICT_MODIFIED.shp ELECTORAL_DISTRICT.shp

But, I get number of warnings in ogr2ogr with message, even though the New file is created and :

Warning 1: Value ...... of field structural of feature ..... not successfully written. Possibly due to too larger number with respect to field width

Afterward, in both cases, I confirm the EPSG in QGIS. It shows the EPSG:4326, WGS84.

But, when I try to load this converted shapefile to local DB in SQL Server, using the below-mentioned command:

ogr2ogr -f MSSQLSpatial "MSSQL:server=RBANSAL\SQLEXPRESS;database=MySpatial;trusted_connection=yes" C:\Users\rbansal\Desktop\ShapeFile\FED_CA_2_2_ENG_Modified.shp -nln "Spatial_Data" -progress

I get the following error.

ERROR 1: INSERT command for new feature failed. [Microsoft][ODBC SQL Server Driver][DBNETLIB]ConnectionWrite (send()).
ERROR 1: Unable to write feature < feature no > from layer < layer name >.
ERROR 1: Terminating translation prematurely after failed translation of layer < layer name > (use -skipfailures to skip errors)

I need help with 2 questions,

Q1: Am I converting the EPSG:5320 to 4326 correctly?

Q2(a) If not, is there any other way?

Q2(b) If I am converting it correctly then why am I getting the error while loading the shapefile to SQL Server?

I tried as per @DPSSpatial's solution:

ogr2ogr -overwrite -nln "Spatial_Data" -f MSSQLSpatial "MSSQL:server=RBANSAL\SQLEXPRESS;database=MySpatial;trusted_connection=yes" "ELECTORAL_DISTRICT.shp" -s_srs EPSG:5320 -t_srs EPSG:4326 -lco geom_name=shape -lco UPLOAD_GEOM_FORMAT=wkt

But, again, got the same error:

ERROR 1: INSERT command for new feature failed. [Microsoft][ODBC SQL Server Driver][DBNETLIB]ConnectionWrite (send()).
ERROR 1: Unable to write feature < feature no > from layer < layer name >.
ERROR 1: Terminating translation prematurely after failed translation of layer < layer name > (use -skipfailures to skip errors)


texts:
Entry layers:
polygons - no attributes (columns/values), points - needed attributes (columns/values)

How to get this result?: polygon layer with attribute (column/value) of most often (most frequent - MODUS) of values from points values within polygons

Probably best for me with SQL select in Mapinfo but I can convert layers and do it in Qgis or ArcGIS if you know the simple way in these SW.

Additionally: How to get count of points within polygons and add it as column (attribute) to polygon layer in Mapinfo? - solved - this not help to main question

Example of points within polygons - polygons represented by outlines (a lot of overlapping - it is import from "nonGIS" software through dxf - nobody cared about topology etc.):

Here is example how to get maximum/minimum/mean... there is a simple solution in ArcMap - right click on layer - Joins and Relates / Join... / Join data from another layer based on spatial location (as described here: In ArcGIS how can I append the number of points that are located inside each individual block groups of the entire US?) - but I need MODUS.


texts:
I am using QGIS3.4 (built-in Python 3) actions.addAction and trying to open a local photo, with the below code I am able to store the photo path in the feature attribute and created an action, I have also ticked the box "auto open form" in identify results tool

QgsProject.instance().addMapLayer(vectorLyr)
iface.setActiveLayer(vectorLyr)
activeLyr=iface.activeLayer()
actions=activeLyr.actions()
actions.addAction(QgsAction.OpenUrl,"Photos",'[% "photo" %]')

I verified the action has been created successfully:

The problem is when I click on the feature using identifier tool, the popped up window is not the photo, instead showing the photo below (I have also tried to add "file:///" in front of the path attribute, and still got same issue)

Then I also tested other type of actions (this is related to my second question below) such as GenericPython to see if it works using below line, again when I clicked the feature, nothing happened:

actions.addAction(QgsAction.GenericPython,"Test Action","print(1)")

I am not sure what I missed here and didn't quite understand the usage of action scopes and my actual goal is to define the action to be able to call another python script (.py) how do I modify the code to make it work?


texts: Running Python script that includes ArcPy functions?

texts:
I have to export about 250 layers in one mxd File to maps in ArcGIS 10.

I found a script on this site at Exporting each layer in map to separate image using ArcPy?:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, '')[0] # assuming there is only 1 df you're interested in

#Read input parameters from script tool
allLayers = arcpy.GetParameterAsText(0) # Must be a multivalue in script tool params
lyrList = allLayers.split(";")

PNGPath = arcpy.GetParameterAsText(1)

#Turn of all lyrs in list
for lyr in arcpy.mapping.ListLayers(mxd, '', df):
for layer in lyrList:
    if lyr.name == layer:
        lyr.visible = False
arcpy.RefreshActiveView()

# Loop through each layer, turn it on and export map as PNG
for lyr in arcpy.mapping.ListLayers(mxd, '', df):
for layer in lyrList:
    if lyr.name == layer:
        lyr.visible = True
        arcpy.RefreshActiveView() # May want to test without this -- ArcMap might export correctly without need for refresh active view
        arcpy.mapping.ExportToPNG(mxd, PNGPath+"\\" + lyr.name + ".png")
        lyr.visible = False

I never before handled Python.

I opened Geoprocessing -> Python in ArcGIS and copied the code wihtout changes. So far as I understand it nothing happened... Which part of the code do I have to change or add to make it run on my computer?


texts: ArcGIS project raster mosaic bends the raster but not the footprint when projecting

texts:
I have a mosaic whose rasters are using a custom coordinate system:

US_IFR
Authority: Custom

Projection: Lambert_Conformal_Conic
false_easting: 0.0
false_northing: 0.0
central_meridian: -95.0
standard_parallel_1: 33.0
standard_parallel_2: 45.0
latitude_of_origin: 39.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_North_American_1983
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_North_American_1983
  Spheroid: GRS_1980
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314140356
    Inverse Flattening: 298.257222101

After adding the rasters, I modify their footprints to trim off unnecessary parts using a feature class to replace the original footprints.  The feature class was created in the same custom coordinate system.  At this point, the mosaic looks great, all the images are properly stitched together and all the extraneous information is removed.

However, once I project to web mercator:

WGS_1984_Web_Mercator_Auxiliary_Sphere
WKID: 3857 Authority: EPSG

Projection: Mercator_Auxiliary_Sphere
False_Easting: 0.0
False_Northing: 0.0
Central_Meridian: 0.0
Standard_Parallel_1: 0.0
Auxiliary_Sphere_Type: 0.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
  Spheroid: WGS_1984
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314245179
    Inverse Flattening: 298.257223563

in preparation to publish as a map service, the rasters get bent, but the footprints remain straight lines.  This means that some of the desired content is stripped off and some of the undesired content is let in.

Any idea why the footprint did not bend in the same way the raster did?

US IFR:

Web Mercator:


texts:
I have used openlayers along with geoserver. My Map is having LCC projection by default, but SRS is EPSG:4042.

The problem is that the scaleline shows NaN in meter as well as feet. How can I correct it?

var bounds = new OpenLayers.Bounds(
    66.76242541918388, 6.756367055683252,
    98.64805108111362, 37.558217731396034
);

var map = new OpenLayers.Map('map',
    {
        maxExtent: bounds,
        maxResolution: 0.10,
        projection: "EPSG:4042",
    center: new OpenLayers.LonLat(79, 21)
     }
);

EDIT 1 :
This happens only when I give a projection for OpenLayers.Map. Otherwise, the scale is correct.

EDIT 2: This worked for me:

var bounds = new OpenLayers.Bounds(
    66.76242541918388, 6.756367055683252,
    98.64805108111362, 37.558217731396034
);

var map = new OpenLayers.Map('map',
    {
        maxExtent: bounds,
        maxResolution: 0.10,
        //projection: "EPSG:4042",
    center: new OpenLayers.LonLat(79, 21)
     }
);

i.e. commenting out projection.


texts: How to manually control react-leaflet popups (via props) within a GeoJSON layer?

texts:
I'm using react-leaflet to render a GeoJSON feature collection that has points and linestrings:

I was able to get the click and hover events on the actual features themselves to work just fine. But I'd like to be able to hover on a list item (related to the items on the map) and have the popup open. I've been sorting through the docs, github, and trying different things out. But it seems there's no way to do this. Or that I'll have to render my linestrings and points manually vs. using the <GeoJSON data=

My map works well with the click events:

return (
        <Map
            style={{
                height: '100%',
                width: '100%',
                margin: '0 auto'
            }}
            ref={(el) => {
                this.leafletMap = el;
            }}
            center={position}
            zoom={10}>
            <TileLayer url='https://api.mapbox.com/v4/mapbox.outdoors/{z}/{x}/{y}@2x.png?access_token=pk.eyJ1IjoiYWJlbnMwefwefEiOiJjajJ1bDRtMzcwMDssmMzJydHdvcjF6ODA5In0.xdZi4pmkhj1zb9Krr64CXw' attribution='&copy; Mapbox' />
            <GeoJSON data={locations} onEachFeature={this.onEachFeature} />{' '}
        </Map>
    );

onEachFeature works as it should

onEachFeature = (feature, layer) => {
        console.log('onEachFeature fired: ');
        layer.on({
            mouseover: (e) => this.MouseOverFeature(e, feature),
            mouseout: (e) => this.MouseOutFeature(e, feature)

        });
    };

But I don't see how to call the layer.bindPopup without using onEachFeature. How does one change call these methods using prop values? I'd like to let people hover on a list item and have it toggle the popups.


texts: Feature To Point for very small polygon in ArcGIS Desktop gives ERROR 000072: Cannot process feature?

texts:
In ArcGIS 10.1, I'm trying to run Feature to Point for a polygon shapefile, selecting the box to generate points inside polygons. The process is failing with:

  ERROR 000072: Cannot process feature with OID 98720

The documentation suggests I run Repair Geometry, which I'd already done. I've checked the offending feature; it's very small, but has vertices. I can't see anything wrong with it.

Any idea what might be causing the error?


texts: Assigning Raster Generated in Function as Global Variable with ArcPy?

texts:
I am trying to write a script to be as robust as possible and am trying to assign a variable to each of my raster outputs (e.g. "winter") so that I may reference it many more times later in the script.

Since "winter" is being generated (and unsuccessfully assigned) within a function, I wanted to make it a global variable so that I can refer to it in the next function.

However, I'm receiving an error saying that my local variable is referenced before assignment and I don't understand why as I believe I am doing this in one of the lines of code below: winter = outTimes.save("Winter")

I am trying to keep my set up of generating these rasters within functons because I am running multiple similar functions and want to keep the variables that create the rasters local. (The raster outputs themselves are what I want to be global).

import arcpy

arcpy.CheckOutExtension("Spatial")
arcpy.env.workspace = arcpy.env.workspace = """C:\Users\inputData.gdb"""

#---------------------------------------------------------------------------
anthro = arcpy.Raster("Anthro_Features")
conifer = arcpy.Raster("Conifer")
lek = arcpy.Raster("Lek")

def winterHabitat():
    # Multiply Rasters
    outTimes = anthro * conifer

    # Save the output
    global winter
    winter = outTimes.save("Winter")

winterHabitat()

# -------------------------------------------------------
# Pre-Breeding

def breeding():
    # Multiply Rasters
    outTimes = winter * lek

    # Save the output
    global Breeding
    Breeding = outTimes.save("Breeding")

breeding()


texts:
I have GPS NMEA (National Marine Electronics Association) data that I collected being stationary and I can see the displacement of 4 meters just by looking at GPX data uploaded to google maps.

I am trying to see this 4 meters by math. I am using Matlab where I read the data in geodetic degrees and used a function to translate it to ECEF in XYZ. I have done that for both the measurement of GPS and my true location and subtracted both in XYZ.

The result is showing a difference of 312 in X-axis and 60 in Y-axis, instead, I was expecting around 4 for both.

Below is the code I used for Geodetic DD conversion to ECEF XYZ:

  lat = lat/180*pi; %converting to radians
  long = long/180*pi; %converting to radians
  a = 6378137.0; % earth semimajor axis in meters
  f = 1/298.257223563; % reciprocal flattening
  e2 = 2*f -f^2; % eccentricity squared

  chi = sqrt(1-e2*(sin(lat)).^2);
  X = (a./chi +h).*cos(lat).*cos(long);
  Y = (a./chi +h).*cos(lat).*sin(long);
  Z = (a*(1-e2)./chi + h).*sin(lat);

  example of latitude/x-axis allocation:

  GPSDD(1,1) = 50.682201667 >> XYZmeas(1) = 3976168.03008252

  TrueDD(1) = 50.682214 >> XTZtrue(1) = 3975855.46252364

  XYZmeas(1) - XTZtrue(1) = 312.567558876704

Did I miss anything or is there any suggestion what would be my mistake?


texts: Forestry Service Raster ImageServer in QGIS

texts:
I routinely use these Raster Services (http://forestrymaps.alaska.gov/arcgis/rest/services/Rasters) for reference in GoogleEarth, but I'd much rather have access to them in QGIS. I've tried various combinations/that have worked for other services I stream, however this one has me stumped. I tried installing the ArcGIS REST API Connector, however it keeps giving me an error. I've reviewed various posts on this forum and haven't found a solution yet. I am running the stable build 2.14.16-Essen

Combinations I've tried:

Using Layer URLs In ArcGIS REST API Connector:

http://forestrymaps.alaska.gov/arcgis/rest/services/Rasters/AK_63360_Topo/ImageServer/exportImage? (This is the address that my Google Earth Network Link is able to use)
http://forestrymaps.alaska.gov/arcgis/rest/services/Rasters/AK_63360_Topo/ImageServer
Various vector WMS layers from the same site (I get syntax errors when I try to connect these)

Using Other tools:

Standard WM(T)S Add layer tool. (Tried more address links than I can count...)
WCS Add Layer. This at least tries to connect, but I get error windows that vary depending on the address link tried

My understanding of how these services function is mostly gained by trial & error. This is one I really want to use, and I know I can access it through Google Earth, so I must be only limited by my own ignorance.


texts: Two different outputs form QGIS basic statistics

texts:
I have two separate line layers, each with a length attribute column.  When I compute basic statistics on the first layer's column, it works as expected, producing Mean, StdDev, and so on

:

But when I input the second layer's length column, I get a different result:

Both columns are in real number, length 20, precision 2 format. There are no null values.  Why is this happening?


texts: Multiple geometry columns for points, lines, and polygons in SpatiaLite

texts:
I'm rather new to the GIS world and even newer to Spatialite and trying to create a new DB for our application. I would like to create spatial index(es) on the geometry column(s), but our geometry BLOB's have always existed in a single column as WKT. After doing some research, it looks like you need a separate column for points, lines, polygons, multipoint, etc.

I'm assuming this from Recipe #6 in the Spatialite Cookbook which states that you need to create your table first, then use AddGeometryColumn to actually add the column which will contain the geometry. In creating this column, one of the items you specific is the GeometryType.

I would much rather keep a single geometry column that could handle all of the shape types, but if the only way to create spatial indexes is to have separate columns, I'll go that route.

Folks before me on this project created their own RTree tables and index to pull data that was within an MBR. My goal is to get rid of the home-grown stuff and use Spatialite for our queries. We have a number of shape files that contain points, lines, multilines, and polygons and I would rather insert their BLOB (WKT) into a single geometry column vs. 1 for each type. Our main queries are to find all shapes that are within an MBR.


texts: Getting around no data places on Kriging in ArcGIS Desktop?

texts:
I'm working with a natural vulnerability research, using a methodology that interpolate data from a high number of wells (around 600) to create a map of the vulnerability of the area.

Following the methodology, I intepolate the data using Kriging.

Since I have and region with almost no wells in my study area, I would like to know if there is any way get around this situation and don't get a map with the kriging error like on the image.

I'm using ArcGis 10.3.1, the interpolation was made on Raster Tool from Surfer.


texts: Is Cloudmade still updating OSM data downloads?

texts:
I sent this question to Cloudmade but unfortunately they never responded. Hopefully someone here knows something about how they work.

They provide an OSM download for the Americas which I am using in a project. Although their site states that downloads are updated roughly every month it looks like nothing has changed since December 2011.

Although there are many country and continent-specific downloads available for OSM data online it seems like this is the only place I can get something for all of the Americas.

Any input on Cloudmade's update strategy or alternate sources for Americas data (without setting up my own planet-wide PostGIS) would be very much appreciated.


texts: QGIS: Arrow symbol from 4 raster layers (x y dir)

texts:
I have 5m raster layers which represent runoff velocity. They are in 4 separate  layers;
velocity, X-location, Y-location, and Direction.
I want to create an arrow symbol with color according to velocity raster and arrow location and direction from the other 3 layers. How can I do this in Qgis?


texts:
I want to create Analysis report after creating model for generating cut/fill.

Once it will done, based on this cutfill/table, would be display statics in html format using Python.

Is it possible, and can you tell me how it will done?


texts:
I use extjs4.1 with geoext2 in my web app.
I want to load mapPanel to tabPanel via ajax and it fill tabPanel height and width.
If i set mapPanel height and width i can see it in tabPanel but it can't fill tabPanel because i set it height and width.Following code:

mapPanel = Ext.create("GeoExt.panel.Map", {
            renderTo: "mappanel",
            height: 575,
            width: 1124,
            map:map,
            zoom: 11,
            tbar : toolbar
        });
var tabs = Ext.create('Ext.tab.Panel', {
            region: 'center', // a center region is ALWAYS required for border layout
            deferredRender: false,
            activeTab: 0,     // first tab initially active
            items: [{
                title: 'History',
                autoScroll: true
            }, {
                title: 'Center Panel',
                autoScroll: true
            }]
        });

I insert MapPanel to tabPanel when user click on grid column that has icon via following code.

var grid = Ext.create('Ext.grid.Panel', {
        store: store,
        columnLines: true,
        columns: [
            {
                menuDisabled: true,
                sortable: false,
                xtype: 'actioncolumn',
                width: 20,
                items: [{
                    icon   : '<?php echo Yii::app()->baseUrl.'/images/icons/show.png';?>',
                    tooltip: 'Show Map',
                    handler: function(grid, rowIndex, colIndex) {
                        var rec = store.getAt(rowIndex);
                        var jsonData = Ext.encode(store.proxy.reader.jsonData);
                        tabs.remove(tabs.getComponent(2));
                        tabs.insert(2,{
                            title:'Map',
                            layout: 'fit',
                            loader: {
                                scripts: true,
                                autoLoad :true,
                                params:{
                                    history:jsonData,
                                    index:rowIndex
                                },

                                failure : function(){
                                    alert('failed');
                                },
                                url: '<?php echo Yii::app()->createUrl('MapWidget/HistoryMap');?>'

                            }
                        });
                        tabs.setActiveTab(2);
                        tabs.doLayout();
                       }
                }]
            },
...

But when i clean height and width in mapPanel options i can't see it in tabPanel!
How can i load mapPanel to tabPanel to fill it's content?


texts: Automatically clean GPS based line files that have multiple traverses and GPS drift using QGIS 3

texts:
Given a GPS Log with floating lines within ~5m at a time, is there any process that will automatically clean these and smooth these out?

I have tried Smooth, Simplify, Saga Line Smoothing and Cartographic Line Generalization with various options.

None of these get rid of this level of issues unless settings like the following are used

processing.run("native:smoothgeometry", {'INPUT':'Z:/My Drive/Mangoesmapping/Spatial Projects/2019/020_Adhoc_CAFNECRide/Working/2019_Route_Edit_Master/2019_CAFNEC_MGA94z55.shp|layername=2019_CAFNEC_MGA94z55','ITERATIONS':1,'OFFSET':0.5,'MAX_ANGLE':90,'OUTPUT':'TEMPORARY_OUTPUT'})

But this gets rid of the corners as well and makes the line to smooth.

Ideally I just want all the nodes between the measure tool start/end to be removed automatically and the remaining nodes to be about where the measure line starts/stops.

UPDATE 1

Based on suggestion from @babel I tried Multipart to Singlepart -- > $length --> Select < 0.5 and this does a good job of showing the problem areas.

If I now delete these I get gaps in the data and also some slightly longer slivers left

Dissolving them retains the gaps -is there a way to remove the gaps programmatically.

I've uploaded a dataset to test with at https://drive.google.com/drive/folders/1eaXpIZvUfMhzEtnGmHe9AMgH2ClQBfiX?usp=sharing


texts:
I am using Kansas Cropland raster data from NASS with tabulate area to find out acreage of different crop types.  The python I used works fine for 2006-2009 (56 meter cell size), except 2010 & 2011 when they switched the cell size to 30 meter.  Every time I ran it would generate empty attribute table (zero values for all field).  I checked raster dataset properties to compared 2010 with previously successfully ran rasters, the only different is the cell size and not allowing null value after 2010; attribute tables setup are identical.

The python I used to ran tabulate area, as shown following:

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "H:/Projects/Ogallala_2012/Docs/Tables"

# Set local variables
inZoneData = "H:/Projects/Ogallala_2012/Geodatabase/Imagery.gdb/ks_counties_30m"
zoneField = "FIPSSTCO"
inClassData = "H:/Projects/Ogallala_2012/Geodatabase/Imagery.gdb/NASS_cropland_2010"
classField = "Class_Name"
outTable = "H:/Projects/Ogallala_2012/Docs/Tables/cropland2010"
processingCellSize = 30

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute TabulateArea
TabulateArea(inZoneData, zoneField, inClassData, classField, outTable,
             processingCellSize)

I've tried to use both raster and feature data to define the zone area,but attribute table still came back empty.  However, when I ran a zonal statistics, it came back with values that seems right.  I also ran tabulate area in both model builder and toolbox it would just crash.

Attached is the empty attribute table it generated.

Thanks for helping!


texts: Replacing Null attribute with the value of the nearest feature

texts: I have several polygons with "Null" values that I would like to give them the "Name' of their nearest polygonal feature. As you see in the photo, the light green polygon is actually a part of the purple polygon, but it is separated because it has null name and not the same name as the parent polygon. Any idea how can I assign the same name to both polygon in order to end up with one harmonised feature using FME, Please?


texts: How to retrieve the attributes of the clicked line or polygon feature in GeoServer?

texts:
In GeoServer, I'd like to retrieve the attributes of a feature that is clicked by the user. When the feature is a point, I can simply use the bounding box parameter of a WFS GetFeature method to get the point.

I do not know how to do when the feature is a line or polygon. Any ideas? An example would be perfect.


texts:
I have a shapefile created by a corn planter. It creates a point for each planting row every second. So a 16-row planter creates 16 points as it goes along the field. As a consequence, a small 30-acre field has over 400,000 points and it's a nightmare to load, not to mention to start any analysis.
The points are created systematically in increasing order. I would like to extract points every Nth observations. For example, to get the center of the planter pass, I would extract the FID=4, and then every 8 observations after that. Any suggestions to do this?

I'm using ArcGIS 10.2 for Desktop.


texts: Converting point XY to page units XY using ArcPy

texts:
I am using ArcPy to export PDFs from my ArcGIS Server application. I want to move a TextElement in my layout view to the centroid of a selected feature. I can get the centroid of the feature, but that's in my projected coordinate system (PCS). I need to convert that point into inches in my layout view.
elements = arcpy.mapping.ListLayoutElements(mapDoc, "TEXT_ELEMENT")
    for elem in elements:
        if elem.name == "SiteText":
            elem.elementPositionX = center.X
            elem.elementPositionY = center.Y

center.X and center.Y are in a projected coordinate system (PCS), not inches. I know it can be done with ArcObjects, as I've done it in desktop addins using C#. I need to do it from ArcPy and I don't see how. Any ideas?


texts: Vector layer parameter in QGIS script takes wrong layer

texts:
I have 6 postgis layers in my QGIS project and all of them have the same data source but different joined data.

When accessing these layers from a python script with a vector layer parameter (like this ##layer=vector) it will always take the same layer (the first one in alphabetical order), no matter which one I select in the drop down selector.

   ##layer=vector
   import processing
   lyrName = processing.getObject(layer)
   progress.setInfo(lyrName.name())

I don't know what I am doing wrong. Any idea?


texts:
Is possible to set up web access to GIS data such that a user could login and see the collected data and be allowed to suggest/draw changes but not be able to edit the data files?

My idea is that community members could see data that has been collected through land use interviews, then be able to suggest changes to the project maps through a web based editor. A GIS tech would then review the suggestion to see if the change would be appropriate.

One issue we have is that our data is considered highly proprietary and confidential. So security is very important.


texts: Openlayers Vector Layers eventListeners: how to specify different data without defining multiple layers

texts:
I have a map. It has two origins, and two destinations. The code for one origin is here:

vectors_origin = new OpenLayers.Layer.Vector("Origin", {
            eventListeners:{
                    'featureselected':function(evt){
                    var feature = evt.feature;
                    var popup = new OpenLayers.Popup.FramedCloud("popup",
                    OpenLayers.LonLat.fromString(feature.geometry.toShortString()),
                    null,
                    json.orig1_iata_code.iata+"
",
                    null,
                    true,
                    null
                    );
                    popup.autoSize = true;
                    popup.maxSize = new OpenLayers.Size(400,800);
                    popup.fixedRelativePosition = true;
                    feature.popup = popup;
                    map.addPopup(popup);
            },
            'featureunselected':function(evt){
            var feature = evt.feature;
            map.removePopup(feature.popup);
            feature.popup.destroy();
            feature.popup = null;
            }
    }
    });

Note this line:

json.orig1_iata_code.iata+"
",

This takes care of one origin. Now since all 4 points (2 origins, 2 destinations), this HTML is going to differ a bit. There will be an

orig1_iata_code, orig2_iata_code, dest1_iata_code, and dest2_iata_code

I really don't want to have to define 4 separate vector layers. Sure it will work, but that seems pretty tacky (correct me if I'm wrong).

So where that html is called above, json.orig1_iata_code.iata+"
",, is there any way to specify 'multiple' optional parameters? Multiple I've done before, but all of them will be called every time.


texts: How to change Geoserver Mapfish print option from pdf to gif


texts: Arcpy 'Append_management' much slower than ArcMap Tool?

texts:
Basically I have a script that grabs data from a web service, parses it, stores it in an intermediary local fgdb, then appends that data to the feature class in the enterprise gdb.

The portion of the script that is taking forever is just this call inside of a function:

arcpy.Append_management(source_fc, target_fc, "NO_TEST")

From one featureclass in the fgdb to another, it is lightning quick

From the fgdb to the enterprise geodatabase, it takes 19 minutes. No matter the amount of rows. 10 rows, 100 rows, or 300 rows all take 19 minutes.

I tried the same exact data using arcpy and then using ArcMap. ArcPy took 19 minutes, ArcMap took 1 second. (For 138 records). These were tested on the same machine.

Any reason this would be doing this?

The 'target_fc' in this instance is not registered as versioned.


texts: Area of classified pixels in an image collection

texts:
I am using Monthly History data of the Global Surface Water dataset (Asset id: JRC/GSW1_0/MonthlyHistory).

I have filtered the images using .filterdate and .filterbounds for my region of interest (ROI) and time period. Images are classified into water, no water and no data classes.

Now I am trying to calculate the area under each class in my ROI.

How can I do that?


texts:
I am confused a bit, using addVectorLayer() in QGIS plugin.

I want to add and display a vector layer, directly when QGIS opens. Hence I am using addVectorLayer() in ____init____() function of my plugin. The problem I am seeing is, the vector layer is loaded in the legend but is not displayed on the canvas. I, manually have to click on "zoom to layer extent", to display it.

Any hints on how to solve this problem will be of high help to me.

NOTE: If I use either of the below, it works perfect (i.e. the layer is loaded and displayed as well).

qgis.exe D:/sample.shp at the command prompt OR
In QGIS Plugin, If I use addVectorLayer() on some menu item click


texts: How to symbolize a large amount of points with lots of color

texts:
I have an interactive map with a very large amount of different points.  Depending on user selection, various points are shown or not-shown.

These points are divided into several categories, each with potentially 20-30 sub-categories.

My current idea is to create a general hierarchy like the following:

Cat 1 => blues
Cat 2 => greens
Cat 3 => reds
Cat 4 => browns

Each of these color groups would be divided into colors based on some criteria in each category.

I'm not concerned with color blindness issues (too many colors needed to worry about this - plus, there is a popup attached to each point giving far more detail).  I'm not even really concerned too much with differentiation beyond the 4 categories, but I did think it might look better if there was some differentiation.

I've looked at Colorbrewer, but they don't have recommendations for this amount of colors.

So my question is what you all do for color symbology in situations like this?  What are (if any) best practices when you have to deal with a lot of different colors?


texts: QGIS outer glow effect without self-reinforcement

texts:
I need to have an outer glow on the entire coastline of my map. The way QGIS creates glow effects with a gradient from alpha=1 to alpha=0 creates a problem, that whenever the glow overlaps with itself because of the geometry of a shape (e.g. bays or islands near the coast), it will reinforce its own strength at the overlapping location (see image). This creates an uneven distribution of the glow effect.

I noticed that there are a bunch of different blend modes for the outer glow effect, yet none of them seem to prohibit this kind of self-reinforcement.
How do I resolve this issue?


texts: Change the geometry type from LineString to MultiLineString SQL Server

texts:
I have a SQL Server database table to contains two geometry types: LineString; and, MultiLineString.

What is the best way to change the type of all of the LineString geometries to 1 part MultiLineString geometries?

I have tried:

UPDATE Footpaths
SET Geom.STGeometryType() = 'MultiLineString';


texts: Add Excel-style row_number and column_name attributes to a Grid Index in ArcGIS

texts:
After using the Grid Index Feature tool it creates a grid polygon layer. The resulting feature class contains a field called [PageName], which displays like A1, A2, A3...etc in each row starting from left to right. In each cell I notice that it shows that the columns (left to right) are labeled as numbers (1,2,3,4... etc) and the rows (up and down) are labeled as letters (A,B,C...AA,AB... etc).

What I want to do is add attributes so that the Columns are labeled as letters and the Rows as numbers; similar to an Excel sheet. I would prefer to have the results in 3 different fields, one for Columns, one for Rows, and one for the CellID (Columns + Rows). How can I do this?

Edit 1: While trying to create a tool I have this so far for the rows but nothing happens.

import arcpy

# set parameters
fc = arcpy.GetParameterAsText(0)
field = arcpy.GetParameterAsText(1)

with arcpy.da.UpdateCursor(fc, field) as cursor:
    for row in cursor:
        row[0] = row[0]
        cursor.updateRow(row)

Edit 2: Clarified the above question more. My test grid I am working with is 46 Rows by 187 Columns. The field calculator this isn't exactly what I am looking for, but this is the progress I have so far if this helps someone. In a new field (Row) I used this code to generate row numbers.str((!OID!+186)/187)


texts: How can i show spatial data(Point, LineString,Polygon) on bing map using openlayer3?


texts: How does Reduce to List works and how to roll back to an image

texts:
I have used reduceRegion using reducer: ee.Reducer.toList() to create a list of NDVI pixel values of a certain image.
I wish to use these values to input them into a CNN for regression.
When I try to transform the list to an image, I am not sure about the shape that I need to reshape the list into.
For example - I have a region that yields 247 pixels when I reduce the image to a list (scale: 10).
That does not allow me to make it a 14X14 image or 15X15 image so I can use it as input to the CNN (as a np.array(15,15)).
How can I transform back into an array?
Also - what is the method used by the reducer?

Is it going row-by-row of the image or "snake-figure"? That matters to me because if it is row-by-row than at some points in the list, two pixels that are next to each other are actually the end of one row and the start of the next one and thus I cannot use the relationship between them.


texts:
I am using Geonetwork  for metadata storage and search and from what I have found I got to knoew that it does support the ISO 19110 meant for the feature catalog but I want to know where do we add the feature description in the metadata of geonetwork
Please help me with this as I can not find anything in the documentation.


texts: How to cache .mxd without data frame's background?

texts:
I'm trying to cache a map document in ArcGIS in order to upload it later to my ArcGIS Online account. My map document (.mxd) contains several layers - most of them are from the same feature class. The layers differs from one another in their symbology and labeling according to specific extents (I've defined an extent interval for each layer specifically). The feature class projection is defined as GCS_WGS_1984.
In order to do so I've made a 1. Scheme with the Generate Tile Cache Tiling Scheme (Data Management) tool 2. and then tried to create a cache with the Manage Tile Cache (Data Management) tool.

I've encountered two major problems:
1. After creating the cache, I've added it to my map document. The cache had a blocked (not-transparent) background (in the color of the data frame's background).
2. The cached raster that I've added, didn't contained the full extent of the .mxd's layers.

My questions are this:
1. How to make the caching procedure to omit the blocked background?
2. How to make the final cache in the .mxd's full layers extent?

I've tried about 15-25 attempts to do so. They included drawing a specific AREA OF INTEREST, transformations of the projections of the data frame's layer, using a default symbology for the layers (which worked but had a bad symbology), disabling the maplex engine, define a full extent for the data frame and more.


texts: MakeFeatureLayer_management where_clause input by the user

texts:
I want to convert a script into a toolbox and one of the parameters is selecting cities with a population > than 500000 hab. I'm using the function MakeFeatureLayer_management.

How can I do to introduce the value of the {where_clause} ('"POP_MAX" > 500000') using arcpy.GetParameterAsText?

import arcpy
from arcpy import env
env.overwriteOutput = True

wsp = env.workspace = arcpy.GetParameterAsText(0)
wsp_result = arcpy.GetParameterAsText(1)

try:
  arcpy.MakeFeatureLayer_management ("cities" + "layer_cities",'"POP_MAX" > 500000')


texts: Seeking administrative divisions shapefile for Italy?

texts:
I am looking for administrative divisions' shapefile for Italy for QGIS. I found one which is a bit old with some inconsistencies.

Can anyone please suggest a database to find better (or newer) shapefiles?


texts: How to find the maximum-area-rectangle inside a convex polygon?

texts:
In this post we are looking for algorithms / ideas on how to find the maximum-area-rectangle inside a convex polygon.

In the following figure, numbers are the areas of the fitted rectangles. As shown a desired rectangle can vary in each dimension and can be in any angle.

Edit:

  We don't have any clear idea how to deal with the problem mentioned as
  so we asked here. Nevertheless, we guess the maximum-area-rectangle
  may be one of those that has one edge aligned on (not necessarily the
  same length edge, of course) an edge of the polygon.


texts: Features to stay selected by default in QGIS?

texts:
Does anyone know if there's a way to configure QGIS so that when you plot a new feature (i.e. a polygon), it stays selected by default?

If I have to select it separately, and it happens to fall on top of other features, it is difficult to select the right one using the standard 'Select Features' button.


texts: Which NIR bandwidth is suitable to Calculate NDVI?

texts:
Often we encounter multi-spectral and hyper-spectral sensors offering multiple NIR bands. For instance, the WorldView 3 sensor has 8 MS bands, out of them, two bands fall within NIR region.

Near-IR1: 770 - 895 nm

Near-IR2: 860 - 1040 nm

Reference (https://www.satimagingcorp.com/satellite-sensors/worldview-3/)

In case of hyper-spectral imagery, there are number of NIR bands with narrower bandwidth. So, in terms of NIR spectral range, there are many choices which apparently makes it hard to select suitable NIR band to compute spectral indices like NDVI.

What are the conceptual and theoretical basis to select optimal NIR spectral range to compute NDVI?


texts: Instructions for installing Sextante plugin

texts:
As a follow-up to the question I recently asked on PyQGIS, can someone please point me to a document, webpage or thread that provides me with step-by-step instructions as to how to install an external QGIS plugin?

I've read through the  manual's "Loading an external QGIS plugin" section and I have not succeeded in decrypting either of the following two threads:
How to install Openlayers plugin on Qgis 1.8?
or
How to install QGIS 1.8 plugin Georeferencer?

I tried to add the Sextante repository https://github.com/qgis/Quantum-GIS.git in the plugin manager, but it remains empty under the plugins tab.  I'm obviously missing some fundamental piece here.  I also cannot find my plugins directory (C:\Program Files\Quantum GIS Lisboa... ????).

Any help is much appreciated,

thanks, mike


texts: Can i set a specific transformation, for my Feature Dataset?

texts:
Using ArcMap, i created a 'Feature Dataset (FD)' and set it to WGS84, F23s.

If i export a shapefile in SAD69 F23s, to this FD, it automatically project the shapefile to WGS, using the "SAD_1969_To_WGS_1984_14" transformation (WKID 1877).

I want to make the same operation, but instead of using this transformation, using the "SAD_1969_To_WGS_1984_13" (WKID 1876).

Is it possible?


texts: R/sf - Transform a geometry column made of POINTs and MULTIPOINTs into xy-coordinates columns (x1, y1, x2, y2, ...)

texts:
After intersecting a set of routes (LINESTRING) with a perimeter (say the border of a country) using st_intersection from sf package, I get a list of points and multipoints: there is one point when the route crossed the border once, and there are many points (multipoints) when the route crossed the border several times. My goal would be then to transform this output into columns of xy-coordinates, with NA when necessary.
Example of input:
library(sf)
g = st_sfc(st_multipoint(matrix(1:6,,2)),
           st_point(c(1, 2)),
           st_multipoint(matrix(11:18,,2)))
df = data.frame(a = 1:3, b = 5:7)
st_geometry(df) = g

With df being:
  a b                       geometry
1 1 5 MULTIPOINT ((1 4), (2 5), (...
2 2 6                    POINT (1 2)
3 3 7 MULTIPOINT ((11 15), (12 16...

The desired output would look like:
  a b x1 y1 x2 y2 ...
1 1 5  1  4  2  5 ...
2 2 6  1  2 NA NA ...
3 3 7 11 15 12 16 ...


texts: Merging some polygons when creating Virtual Layer

texts:
I have a polygon layer 'lot' containing an attribute "lot_number". Then, I want to create a Virtual Layer that merges some of the lots while keeping others separated, as the diagram shows below.
Specifically, in the virtual layer, I want to merge the lots with "lot_number" 1, 2 and 4 into one polygon, and ones with "lot_number" 5 and 6 into another polygon. How could I do that?


texts: Accessing PostgreSQL/PostGIS tables as inputs into QGIS (PyQGIS) standalone processing scripts

texts:
With the help of Stack Exchange I have been able to set up a PyCharm environment in which I am able to run QGIS geo-processing queries, so away from the Python console in QGIS itself. If I point to the directories in which my shapefiles or geopackages are located I can fairly easily access those and generate outputs to an output directory of my choice.
However, I would like to be able to access geospatial tables directly from a PostgreSQL/PostGIS database and use the QGIS based processing libraries to run analyses involving these datasets. I have been trying to find a way to achieve this using

the psycopg2 python library ..and/or
a PostGisDBConnector function that some posts suggest used to work prior to QGIS 3.10 and now seems not to work

In fairness, I am experimenting and haven't found any solution yet. Effectively I want to find a way of essentially pointing to, for example, a line feature layer in the PostgreSQL DB and to a polygon layer such that I get a URI for both that I can then use in the parameters of a processing function such as qgis:intersection.
At the moment I cannot find a method or appropriate syntax to achieve this. Has anyone out there tried to do similar with QGIS 3.14 and Python 3?


texts: Understanding the MUSYM field on SSURGO dataset?

texts:
What does the MUSYM field on the SSURGO data set mean? I have looked on the nrcs website have only found that it identifies soil.

In the data set that I am looking at there are values like Ub and DaB. These do not look like some randomly generated identifiers. They seem like they must have some acronym meaning to them.

Where can I find a look up table for these values or is there a system for deciphering them?


texts:
I have a polyline layer. I want to split this polyline layer at CROSSING POINT.
With two steps I want to do this.

I show the user the CROSSING POINT of polyline layer.
With this CROSSING POINT, I want to split the polyline layer.

I am on Arcgis 10 with .NET/arcpy.
Second step can be done with  SPLIT LINE AT POINT tool. But how can I identify the crossing point?

I have attached image to clear my problem.


texts: Pyshp 1.2 example from documentation fails

texts:
I must be doing something really stupid because I am trying to run an example straight from the PyShp (shapefile.py) 1.2 documentation with the only change being the input shapefile and it is failing. Note that it is 32bit Python 2.6 on Server 2008 R2 64bit.

The code:

import shapefile
e = shapefile.Editor(shapefile="F:/data/IntermediateWorkspace/5day_pgn.shp")
e.poly(parts=[[[5.1,5],[9.9,5],[9.9,1],[7.5,3],[5.1,1]]])
e.record("Appended","Polygon")
e.save('F:/data/IntermediateWorkspace/5day_pgn')

The error:

  Traceback (most recent call last):
    File "F:\scripts\shapetest.py", line 2, in
      e = shapefile.Editor(shapefile="F:/data/IntermediateWorkspace/5day_pgn.shp")
    File "C:\Python26\ArcGIS10.0\lib\shapefile.py", line 1048, in init
      self.records = r.records()
    File "C:\Python26\ArcGIS10.0\lib\shapefile.py", line 520, in records
      self.__dbfHeader()
    File "C:\Python26\ArcGIS10.0\lib\shapefile.py", line 443, in __dbfHeader
      fieldDesc = list(unpack("<11sc4xBB14x", dbf.read(32)))
  error: unpack requires a string argument of length 32

Thanks for any help!!


texts:
I would like to dissolve a feature class by a field. When I run the code, I receive an error that says, "invalid topology [out of memory]" I have 16GB of RAM and terabytes available on the hard drive.

Here is my code:

#Import system modules
import arcpy
from arcpy import env

# Set environment settings
env.workspace = "Q:/LandTrendr_04/lt_ancillary/landtrendr_file_geodatabase/templates_testing/LandTrendr_GRPO_python.gdb/"

# Set local variables
inFeatures = "raster2polyg"
outFeatureClass = "polydissolve"
dissolveFields = ["UNIQUE"]

# Execute Dissolve using UNIQUE as Dissolve Fields
arcpy.Dissolve_management(inFeatures, outFeatureClass,dissolveFields,"","MULTI_PART", "DISSOLVE_LINES")

I would like to dissolve on the 'UNIQUE' field, no stats, create multipart polys, and dissolve lines.

I've tried to run this same tool through the toolbox and it works just fine. I would like to have this in a python script to help with my workflow.

Here is the output from the cmd prompt when I run the script.

Thanks

-al


texts: retrieveTravelModes returning travel modes not present in my network dataset

texts:
I am observing a very strange behavior. I have published a network service on version 10.61.
When I click on

.../NAServer/Route/retrieveTravelModes?f=json

I get following travel modes which are not defined in my network data set.

Walking Time
Rural Driving Distance
Driving Time
Driving Distance
Walking Distance
Rural Driving Time
Trucking Time

When I try to get travel mode on the same network data using

arcpy.na.GetTravelModes

I get following travel modes which I defined in my network data set.

Avoid Planned Roads

I have tried publishing same network dataset twice on my server but I am observing the same behavior.

What I am missing here ?


texts: ArcObjects.NET, append with automatic field mapping not working

texts:
I got a strange problem with ArcGIS 10.4 using ArcObjects.NET. By using the code below, if I want to use the Append tool, I have to send it the source and the destination. Both of them are of Object type. So, I can send paramters of type:

IFeatureLayer

IFeatureClass

ITable

String (being the file path)

If I use the files paths, the field mapping is done automatically. But strangely, if I use an IFeatureLayer, IFeatureClass or ITable instance, the fields are not mapped. Even for the ITable, I get all the rows with empty values...

Is it normal? Of course, I can get a workaround by using the files path but when I use a memory FeatureClass, I can only use the IFeatureClass object because there is no file, so, no file path.

Dim gp As New ESRI.ArcGIS.Geoprocessor.Geoprocessor

Dim lAppend As New ESRI.ArcGIS.DataManagementTools.Append(source, destination)

lAppend.schema_type = "NO_TEST"
gp.Execute(lAppend, Nothing)


texts: What is returned when Select event fires in Openlayers?

texts:
I'm quite new to Openlayers and jscript in general. Using Openlayers 3.9.0.
Can anyone hint where to look up the structure of the 'e' object which is returned when Select event fires?

select.on('select',function(e){
  //...
});

For example I come across such things as e.target.getFeatures(). What else is there in the 'e' object?