texts: Installing NVS Vector Stream Tool with ArcGIS Desktop version later than 10.0?

texts:
I can't seem to get a response form NVS about their Vector Stream Tool add-on, but they did send me a fixed-up version that I have now been able to install.

However, I cannot find any associated toolbox.

I have tried System Toolbox, My Toolbox, looking in the installation directory, and searching in my default geodatabase, but I can't find a way to add this tool (or toolbox) to my GIS.

Has anyone else been able to successfully install the Vector Stream Tool into ArcGIS 10.1 SP1 for Desktop (ArcMap), and might anyone have suggestions?


texts: OpenLayers 3 - ZoomSlider outside the map

texts:
I am building a map with OpenLayers and I want to display a custom ZoomSlider outside of the map.

I know that displaying controls outside of the map goes like this:

JavaScript:

var ctrl= new ol.control.Zoom({
        target: document.getElementById('ctrl')
    });
map.addControl(ctrl);

Html:

But it doesn't work with the ZoomSlider Control...

Any ideas?

Here is a Fiddle which shows the problem


texts: Integrating Canadian WMS with ArcGIS Server?

texts:
We are trying to load a Canadian WMS to our ArcGIS Server and I am told the URL should have "arcgis" or "rest" in it for it to work.

This WMS-Toporama link works find in ArcMap with "Add WMS Server" tool but our web site person cannot seem to load it with our ArcGIS Server. I am not familiar enough with the web site side of things to give a good answer to this problem.


texts: How to fix "Parallel Arrays Must Have Same Length" error for chart

texts:
I'm using a transect with two types of data to collect temperature and aerosol index on Google Earth Engine. I want to print a line chart using multiple dates, but I keep getting two error listings:
Array (Error) Dictionary.toArray: Parallel arrays must have the same length.
Error generating chart: Dictionary.toArray: Parallel arrays must have the same length.
This is the section of the code which is probably the problem.

var distance = startingPoint.distance(500000);
var image =   distance.addBands(summer2018).addBands(winter2018).addBands(summer2019).addBands(winter2019).addBands(aerosum2018).addBands(aerowin2018).addBands(aerosum2019).addBands(aerowin2019);
// Extract band values along the transect line.
var array = image.reduceRegion(ee.Reducer.toList(), transect, 1000)
                 .toArray(image.bandNames());
                 print (array, 'array');

// Sort points along the transect by their distance from the starting point.
var distances = array.slice(0,0,1);
array = array.sort(distances);

// Create arrays for charting.
var aerosoltemp = array.slice(0,1);  // For the Y axis.
// Project distance slice to create a 1-D array for x-axis values.
var distance = array.slice(0,0,1).project([1]);

// Generate and style the chart.
var chart = ui.Chart.array.values(aerosoltemp, 1, distance)
    .setChartType('LineChart')
    .setSeriesNames(['Summer 2018','Winter 2018','Summer 2019','Winter 2019', 'Aero Summer 2018', 'Aero Winter 2018', 'Aero Summer 2019', 'Aero Winter 2019'])
    .setOptions({
      title: 'Air Pollution and Temperatures Along Indonesia Transect',
      vAxes: {
        0: {
          title: 'Average Seasonal Temperature (Fahrenheit)'
        },
        1: {
          title: 'Aerosol Index',
          baselineColor: 'transparent'
        }
      },
      hAxis: {
        title: 'Distance from Jambi (m)'
      },
      interpolateNulls: true,
      pointSize: 0,
      lineWidth: 1,
      // My chart has two Y axes: one for temperature and one for air pollution.
      // The Visualization API allows us to assign each series to a specific
      // Y axis, which we do here:
      series: {
        0: {targetAxisIndex: 1},
        1: {targetAxisIndex: 0},
        2: {targetAxisIndex: 0}
               }
      });

print(chart);
Map.setCenter(103.61, -1.58);
Map.addLayer(aerosol, {min: -1, max: 2.0},'pollution');
Map.addLayer(transect, {color: 'blue'},'transect');


texts: R - Crop a pixel image produced by spatstat with a boundary shape

texts:
i am using the very good spatial statistic package spatstat to produce probability maps of spatial point patterns.
Here is a simple working result from the examples to try out.

  require(spatstat)
  data(lansing)
  z <- relrisk(lansing)
  plot(z, main="Lansing Woods")

My question is rather simple (or maybe not). I produced a similar map as shown above and i want to crop it in R with a boundary shapefile. The spatstat package uses its own system (of windows and pixel images) so i am rather clueless. I know that i can format my shapefile to a so called mask with this command as.mask(as.owin(shapefile)), but i didn't get beyond this point.
Anyone knows how i can crop the produced image with this mask?

EDIT:
Ahh, sry. I already found a solution around it, which works for me. I simply have to make a subset out of the spatial points before converting it to a ppp object.

library(rgdal)
st <- readOGR("GIS-Data","boundary_shape")
wp <- spp[-which(is.na(over(spp,st)[,1])),] # spp is my spatial point shape
rm(st)
#Then use this command from the maptools package
library(maptools)
plp <- as.ppp.SpatialPointsDataFrame(wp)

Nevertheless it would still be nice to know if it is possible to crop a spatstat pixel image.


texts: What is an easy to use web GIS stack to use for making fictional maps?

texts:
I am looking for a technology stack that would allow me to have some sort of OSM-driven custom map with as-easy-as-possible editing via a web frontend. I would require this for creating custom fictional maps for roleplaying games, so they wouldn't have any connection to the real world at all.

As for editing capabilities, I'm eyeing openstreetmap.org's current default web editor, iD. But to be honest, I have no idea what software I would need to make it all work together. I also have only limited knowledge when it comes to setting up virtual machines. If a ready-to-go appliance exists somewhere (or something that requires only trivial extra effort), I would be super grateful.

Edit/clarification: by "fictional map" I mean exactly that, a completely made-up map that very likely doesn't re-use a single node from Earth's actual map but instead is made-up from ground zero. I intend to use OSM for this, mostly because of the easily transferable data format, my (limited but still somewhat present) acquaintance with editing and the vast landscape of readily available mapping/rendering tools. I do intend to use most tags as they would be used on OSM, so editor presets are fine as they are for Earth's non-fictional map.

Edit/clarification #2: it should be web-based as we're physically scattered through-out the world. If we only passed an .osm file around, I fear that we could bump into very inconsistent vaguely different versions of said fictional map too easily.


texts: Cannot install new QGIS 3.16 'Hannover' on Mint 20.1 'Ulyssa'

texts:
I have tried to install the latest release version of QGIS (as per today, 3.16) on Mint 20.1 using the same steps from the QGIS page:
sudo apt install gnupg software-properties-common
wget -qO - https://qgis.org/downloads/qgis-2020.gpg.key | sudo gpg --no-default-keyring --keyring gnupg-ring:/etc/apt/trusted.gpg.d/qgis-archive.gpg --import

sudo chmod a+r /etc/apt/trusted.gpg.d/qgis-archive.gpg
sudo add-apt-repository 'deb https://qgis.org/debian focal main'
sudo apt update
sudo apt install qgis qgis-plugin-grass

However, the system installed QGIS 3.10.4 (an old version), which is the same version available on the software manager.
I have tried to change the repo link to (https://qgis.org/ubuntu) and still did not install the correct version.
Any reason why this is happening?


texts: Is there a bug with the "For" iterator in model builder?

texts:
I have a problem using the "For" iterator in Model Builder (v 10.2.2, 64 bit, Windows 7). I wanted to create a series of rasters using the "Constant Raster" tool. The model works as expected when run from within model builder, but when run as a model tool (both with and without parameters for the For variables) it did the first two iterations, then overwrote the first one for the remaining iterators.

So, with the increment set to 5, the starting timestep as 5 and the final as 20 I would expect to get 3 output layers "Const_5yr", "Const_10yr", "Const_15yr", "Const_20yr". When run from within the model builder, I get what is expected. When run as a stand alone tool, I get "Const_5yr" and "Const_10yr", but "Const_5yr" has the value expected for "Const_20yr" (and from the dialogue in the results window I can see that "Const_5yr" is being overwritten).

To add even more confusion, when I add the iterator model to another model, it works as expected when run as a tool, as well as within model builder.

This is the model builder screen:

I realise this is similar to Does ModelBuilder have Iterate Field Values bug when model run without being in Edit mode? but am not sure if it's due to the specific iterators or a more general problem?  Has anyone else had these issues?


texts: Google Earth Engine - Remove geometry fields from FeatureCollection for export to table

texts:
How do I remove the geometry fields from this feature collection? They are not necessary and give me a messy table in excel.

// Add reducer output to the Features in the collection.
function protected_areas (image, select_country) {var pas = image.reduceRegions({
  collection: select_country,
  reducer: ee.Reducer.frequencyHistogram(),
  scale: 30,
}); return pas.map(function(feature){
  var dict = ee.Dictionary(feature.get('histogram'))
  feature = feature.set(dict)
  return (feature)})
} // appends output from histogram into separate columns

// export function
function export_table (table, description, date) {
  var em = Export.table.toDrive({
    collection: ee.FeatureCollection(table),
    description: region+'_'+description+'_'+date,
    folder: region
});return em}

link to code: https://code.earthengine.google.com/dba0dca5d0395d981ae3e1e97126e49a


texts: QgsMapCanvas object is not callable (Ubuntu with, Eclipse, Python & QGIS 2.0)

texts:
I'm basically borrowing the code for picking up a coordinate from a mouse click from this answer :

How to programatically check for a mouse click in QGIS

However in a stand alone application I do not have access to iface, so am instead calling the class with :

tool = PointTool(self.map_canvas())
self.map_canvas().setMapTool(tool)

Intellisense finds the map_canvas ok and it compiles, but I am getting a run time error:

'QgsMapCanvas' object not callable

Any ideas why?


texts:
I have ArcGIS 10.2.2

I have raster data over a whole state but am looking at just analyzing a specific county. I downloaded county polygons to overlay over the raster data to determine the area I want to study.

I wish to extract the the county information from the whole dataset so I can just work and analyze the information I want.

I tried to use intersect but since I have polygon and raster data (different formats) I couldn't use the command.

Is there another command I could use to extract the data?


texts: Angular distance between Grid and True North given Lat/Long

texts:
As part of a project that uses UTM coordinates I am trying to find a formula to compute the angle (with sign) between the grid north and the true north. I understand this is called meridian convergence.

I am already able to get the correct easting and northing, but I cannot find a reference for a formula that would give me the angle given a single location.

I know that I can derive it mapping multiple points and computing then the angles between the lines that connect them as suggested, for example, here, but if possible I would like to avoid it.

Does such an equation exist?


texts: Band 1: Access window out of range in RasterIO()

texts:
I'm trying to implement a mouse_move event handler to display the elevation based on the mouse location. I'm using the ThinkGeo API to display the raster map (.tif) and gdal C# to get the elevation data from Dted file.
Solve issues that for out of range. But after i play around with the map it suddenly throws :

  ArgumentOutOfRangeException was Unhandled

with message:

  The input double value is out of range.

What this is means?

WORKING CODE(EDIT)

  private void Map1_MouseMove(object sender, MouseEventArgs e)

        {

            Point point = e.MouseDevice.GetPosition(Map1);

            ScreenPointF screenPointF = new ScreenPointF((float)point.X, (float)point.Y);
            PointShape pointShape = ExtentHelper.ToWorldCoordinate(Map1.CurrentExtent, screenPointF, (float)Map1.ActualWidth, (float)Map1.ActualHeight);

            string x = pointShape.X.ToString("f6", CultureInfo.InvariantCulture);
            string y = pointShape.Y.ToString("f6", CultureInfo.InvariantCulture);

            if(proj4Projection != null)
            {
                //PointShape p = new PointShape(605770.08015, 388756.038425);
                PointShape p2 = proj4Projection.ConvertToExternalProjection(pointShape) as PointShape;

                Gdal.AllRegister(); //Register driver(s).

                //string Dted = @"C:\N03.dt2";
                double[] gGeoTrans = new double[6];
                Dataset ds = Gdal.Open(@"C:\N03.dt2", Access.GA_ReadOnly); //Open Dataset

                ds.GetGeoTransform(gGeoTrans);
                OSGeo.GDAL.Band gThisBand = ds.GetRasterBand(1);

                double X = Convert.ToDouble(p2.X);

                double Y = Convert.ToDouble(p2.Y);

                // Your query location is X, Y (both doubles), convert your location into
                // local row, col addressing to specify which cell to read
                int row = (int)((gGeoTrans[3] - Y) / -gGeoTrans[5]);// ULY - xcoord / -N-S CellSize;
                int col = (int)((X - gGeoTrans[0]) / gGeoTrans[1]);  // xcoord - ULX / E-W CellSize

                // check here that the row and col are >= 0 and < gDataset.RasterXSize and
                // gDataset.RasterYSize or you will get a read error.

                // read the value.. reading is done to an array but as we want only one cell
                // it needs to be an array of 1 float.
                float[] buff = new float[1];

                 if ((col <= ds.RasterXSize) && col >= 0)

                 { if ((row <= ds.RasterYSize) && row >= 0)
                    {
                        gThisBand.ReadRaster(col, row, 1, 1, buff, 1, 1, 0, 0); // read the cell value
                     }
                 }
                 else
                 {
                    return;
                 }

                 float ThisCellValue = buff[0];

                TextBoxSatu.Text = "RSO_WEST(KERTAU)" + "(" + (x) + "," + " " + (y) + ")" + " " + "Long: " + DecimalDegreesHelper.GetDegreesMinutesSecondsStringFromDecimalDegree(p2.X) + " Lat: " + DecimalDegreesHelper.GetDegreesMinutesSecondsStringFromDecimalDegree(p2.Y) + " " + "Height :" +ThisCellValue +"m" ;


texts: Transferring data from file geodatabase (.gdb) to cassandra?

texts:
I want to store my spatial data in cassandra.

I want to import data from an existing file geodatabase(.gdb) which has a table containing boundaries of polygons and a table containing point data.

I have individully converted all the tables to csv format and ingested them via fat table to cassandra but I am still worried that the route I am taking is a long one.

Is there a workflow or technique that can directly import .gdb tables to cassandra without converting them to csv?

Storing features in NoSQL database like Cassandra? does not have anything about inserting data.


texts: Finding latitude & longitude values for both ends of a road

texts: I am looking for a way to find latitude and longitude values of both ends of any given road. So for example if there was a site or tool that did this, I would just have to provide a road name (and maybe post/zip code) and it would give me the latitude and longitude values.


texts: FME wont create the data files correctly?

texts:
I try to create the georeferencing files which contains only the one text line like the bounding box coordinates and name (ymin, xmin, ymax, xmax, name) (like: 6635000 3270000 6640000 3275000 18n1a2a).

I create text line using StringConcatenator which connects straight to data writer (Data File). I suppose that it only create file but nothing inside. I am also using AttributeFileWriter in parallel.

This process creates only 100 Data Files. After that it show me error message:

BADNEWS: File C:\Program Files\FME\metafile\MULTI_WRITER.fmf could not be opened (tabrdr.cpp:161)

Anyone who has done similar data files in FME?


texts: Importing EarthVision ASCII to ArcGIS Desktop?

texts:
How do I import a surface map from Petrel in EarthVision ASCII Format to ArcGIS?

I have a file from surface map in Petrel that I want to export to ArcGIS Desktop.

I used the tools 3D ASCII to Feature Class but its shows nothing, just blank.

The tool shows successful processing but when i check the TOC I did not have anything. just blank.

The file type has information about XYZ.

Here is the example of my ASCII file

no header in my file..
its just contains XYZ


texts: Getting blank image from Binary Thresholding function of ArcGIS Pro

texts:
Why do I end up with a blank image every time I run the Binary Thresholding function of ArcGIS Pro on this single band image (Landsat 8, band 4)?
I've run into the same issue before with other raster tools, and changing little things in the symbology like changing the stretch type from percent clip to standard deviation fixed it. But this time I can't seem to resolve the issue and don't know what else to try.
To make sure I was going through the process correctly I used a different but similar image and it worked correctly.


texts: Resample SMOS dataset using cKDTree and Bilinear interpolation

texts:
I have three datasets now. One is SMOS level3 soil moisture dataset which is in EASE grid 2.0 ~25km cylindrical projection. (EPSG: 6933). TRMM precipitation dataset which is in 0.25 degree. SMAP soil moisture dataset which is in EASE grid ~3km cylindrical projection. Now I want to do some data mining work so I want to downscale the resolution of the SMOS and TRMM dataset to be the same as the SMAP dataset. All the datasets are only for Oklahoma state. Say, the SMOS dataset is like a 22 * 43 matrix, now I want to expand it to a 145 * 338 matrix with interpolation with the latitude and longitude being considered.

I am currently trying to use cKDTree in python to find k nearest neighbors for one specific grid. I use the code from this link (the blog only describes the methods but it has an associate notebook with code there). The cKDTree can return four points. Then I can use the latitude and longitude to compute the distance in order to do the bilinear interpolation. Here is the code:

def distance(lat_center, lon_center, lat0, lon0):
    rad_factor = pi/180.0
    lat_center *= rad_factor
    lon_center *= rad_factor
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor

    clat, clon = cos(lat_center), cos(lon_center)
    slat, slon = sin(lat_center), sin(lon_center)
    delX = cos(lat0_rad) * cos(lon0_rad) - clat * clon
    delY = cos(lat0_rad) * sin(lon0_rad) - clat * slon
    delZ = sin(lat0_rad) - slat
    dist_sq = delX ** 2 + delY ** 2 + delZ ** 2
    return dist_sq

def bilinear(lat_center, lon_center, iy, ix, latvar, lonvar, target_var):
    area = []
    for x, y, i in zip(ix, iy, range(len(iy))):
        area.append(distance(lat_center, lon_center, latvar[y], lon_center) * distance(lat_center, lon_center, lat_center, lonvar[x]))
    print area

I do process the data a little bit before construct the KDTree like this:

class Kdtree_fast(object):
def __init__(self, ncfile, latvarname, lonvarname, target_varname):
    self.ncfile = ncfile
    self.latvar = self.ncfile.variables[latvarname]
    self.lonvar = self.ncfile.variables[lonvarname]
    self.target_var = self.ncfile.variables[target_varname]

    rad_factor = pi / 180.0
    self.latvals = self.latvar[:] * rad_factor
    self.lonvals = self.lonvar[:] * rad_factor
    ny = len(self.latvals)
    nx = len(self.lonvals)
    self.latvals = np.repeat(self.latvals, nx, axis=0).reshape(ny, nx)
    self.lonvals = np.tile(self.lonvals, ny).reshape(ny, nx)
    self.shape = self.latvals.shape

    clat, clon = cos(self.latvals), cos(self.lonvals)
    slat, slon = sin(self.latvals), sin(self.lonvals)
    clat_clon = clat * clon
    clat_slon = clat * slon
    triples = list(zip(np.ravel(clat_clon), np.ravel(clat_slon), np.ravel(slat)))
    self.kdt = cKDTree(triples)

def query(self, lat0, lon0):
    rad_factor = pi / 180.0
    lat0_rad = lat0 * rad_factor
    lon0_rad = lon0 * rad_factor
    clat0, clon0 = cos(lat0_rad), cos(lon0_rad)
    slat0, slon0 = sin(lat0_rad), sin(lon0_rad)
    dist_sq_min, minindex_1d = self.kdt.query(
        [clat0 * clon0, clat0 * slon0, slat0], k=4)
    iy_min, ix_min = np.unravel_index(minindex_1d, self.shape)
    return dist_sq_min, iy_min, ix_min

The dataset offers the coordinate information like this:

latitude: 33.04897 33.28237 ... 37.844597 38.092167

longitude: -103.87608 -103.616714 ... -93.24207 -92.98271

I want to implement the bilinear by myself because later I may use some machine learning method to do the interpolation which will still needs to compute the distance between grids correctly.

Now my trouble is the dist_sq_min returned by the cKDTree doesn't match the distance returned by the distance function if I call them manually.

The reason I need to write the distance function is the cKDTree only returns the distance between the target point and the neighbors while the bilinear needs the distance vertically and horizontally to compute the area as the weight.


texts: Convert IEnvelope - ArcMapControl to X,Y screen coordinates

texts:
Does anyone know how to do the following.

Say I have a map and I draw a rectangle over the state of Texas.  ArcMapControl.TrackRectangle returns the lat,lon of this.  However, what I'd like is to actually get the X,Y coords of the screen.  For instance if the map was 500 px wide and 450 tall I'd expect something like

top = 257
bottom = 242
left = 146
right = 167

Anyone know if there is an easy way to convert this?


texts: How Select only specific feature from selected features of more than one layer

texts:
I selected features using select feature tool in a map. I want to fetch only specified features among all selected features of different layers. For example, I want to fetch only roads among all selected features(roads, river, buildings..etc). How to do it in .net?
Please help me out.

Thanks


texts: Adding field to PostGIS Layer with PyQGIS

texts:
I'm having trouble adding a field to a PostGIS table using PyQGIS. Here's my code:

from PyQt4.QtCore import QVariant
from qgis.core import *

inLayer = iface.activeLayer()
provider = inLayer.dataProvider()

# add field
field = QgsField('test_field', QVariant.String)
inLayer.startEditing()
provider.addAttributes([field])
inLayer.updateFields()

# commit
inLayer.commitChanges()

The commit fails and the PostGIS log shows an error:

  2015-10-27T08:13:28    1   Erroneous query: ALTER TABLE "received"."env" ADD COLUMN "test_field"  returned 7 [ERROR:  syntax error at end of input

  LINE 1: ALTER TABLE "received"."env" ADD COLUMN "test_field"

Seems the call to PostGIS fails to include the column type. Is this a bug, or am I going about this the wrong way?


texts: QGIS boundaries of many collected polygons

texts:
QGIS 3.4 LTR
I have a layer with many polygons.
I want 2 symbols for them in this layer :

The first is one color for all (ok...).
the second is the boundary of the polygons grouped by an attribute, using the generator of geometry
So, Here's the start :

I would like

I tried
exterior_ring(collect($geometry, group_by:='nom_dept')) and    boundary(collect($geometry, group_by:='nom_dept'))
But it does not work. Any idea ?


texts: Topographic tile cache needed to replace OSM street map data in an app

texts:
I have developed and marketed an application that successfully uses a pre loaded tile cache of Open Street Map .png files for a small city area about 20km x20km to support the presentation of tourist information on the web, iPad and iPhone. It is written in HTML5 and uses MoSync to create the wrapper to deliver it on the mobile devices.

This approach works well for cities but I would like to extend the app to work in countryside areas (in the UK initially) where if I could access it, Ordnance Survey mapping 25000,50000 and 250000 would deliver a much better user experience.

My specific questions are:

How do I get hold of topographic tiles?
Can I get tiles that have the same slippy map structure as OSM.
Can I use Ordnance Survey as a source for the UK and if so is it a one off license fee?
What zoom levels are available? I expect 25000, 50000 and 250000. Is 100000 available?
What formula do I use to identify the specific tile from a lat and long pair?
Will all the tiles at a specified zoom level be the same delta lat and delta long? If not how do I tile them?


texts: Python Function for Altering shapefile's attribute datatype

texts:
I want to merge a few hundred shapefiles depending on geometry type, which are in several directories and subdirectories. i have written a python script and i am using mainly arcpy for the routines. But the process cannot be completed as there is a problem on the datatype of a common attribute in all shapefiles. In most of shapefiles the datatype of this attribute is string and in some other is double. I found the shapefiles with attribute datatype double, but these are a few, and cannot be done one by one. So i am looking for a function that changes the datatype of this attribute, filtering the current datatype.
Thanks in advance.


texts: DSAS APP doesnt recognize my baseline shapefile says its not projected or the projection is not in meters

texts:
I have created a personal geodatabase which includes a baseline feature layer and a shoreline feature layer both in WGS84 coordinate system. DSAS casts transects in order to calculate shore erosion.

When I try to cast transects an error pops up saying:

Baseline layer is either not projected or the projection does not use metres as linear unit. Edit projection baseline


texts: Identifying and removing bad lines in a shapefile

texts:
I downloaded some shapefiles of Houston from the Open Maps data using this tool. I am plotting the roads and I have a weird thing happening. Most of the lines are correct, but there are clearly several erroneous lines that are drawing long streaks across the map.

How can I identify these lines and remove them from the shapefile?
Here is the code I used to plot the thing.
import geopandas
import matplotlib.pyplot as plt

fname = 'planet_-96.09,29.269_-94.616,30.18-shp/shape/roads.shp'

fig = plt.figure(dpi=300)
ax = plt.subplot(1, 1, 1)

lat_0 = -95.579
lat_1 = -94.99
lon_0 = 29.577
lon_1 = 30.19

ax.set_xlim([lat_0, lat_1])
ax.set_ylim([lon_0, lon_1])

roads = geopandas.read_file(fname)

roads.plot(ax=ax, lw=.05)

plt.savefig('houston.png')


texts: Write attribute column names (fields) in a PyQGIS list based on their attribute value

texts:
I have a tricky question: for a QGIS script, I need to filter a list of attribute column names (fields) of a vector based on the attribute value of each column.
My Vector is structured like this:

I create the list of field names with:
layer = QgsVectorLayer(Input_Layer, "ogr")
fields = layer.pendingFields()
field_names = [field.name() for field in fields]

So field_names gives me a list of all the attribute field names:
[u'NAME_A', u'NAME_B', u'NAME_C']

I now need to filer the list, so attribute fields are only listed in field_names if the attribute value (there is only one attribute) is not NULL. So the result should only list:
[u'NAME_A', u'NAME_C']

I can't figure it out. Running QGIS 2.18.18 on Windows 10


texts: Google Earth Engine Filter seasons independent from year

texts:
Is something like this possible? and how ?

Image.Collection.filterDate('20xx-04-01','20xx-09-31')


texts: Converting decimal degrees to state plane in ArcMap?

texts:
I'm working in Arcmap 10.1 and I have a point feature class that I already have xy values in decimal degrees. The projection is state plane. The feature class also has two fields for the state plane that I need converted from the decimal degrees.

How can I do this?


texts: How can I set Changed flag in ArcMap using ArcObjects.NET?

texts:
In my task I changes representation geometry in selected features.
When I edits something in ArcMap manually, question dialog about saving changes causes:

Editor - Start editing
Edit geometry
Editor - Stop editing
Question Do you want to save your edits?

But when I edits programmatically using ArcObjects.NET, editing stopping doesn't causes question, and all changes losts if they didn't save before stopping.

Where and what flag must I set for question about saving will be caused?

My code for editing:

for (int j = 0; j < map.LayerCount; j++)
{
    IFeatureLayer mapLayer = (IFeatureLayer)map.Layer[j];
    IGeoFeatureLayer geoLayer = (IGeoFeatureLayer)mapLayer;
    IFeatureClass mapFeatureClass = mapLayer.FeatureClass;
    ICursor mapCursor = (ICursor)mapFeatureClass.Search(filter, false);
    IFeature mapFeature = (mapCursor as IFeatureCursor).NextFeature();
    if (mapFeature != null && geoLayer.Renderer is IRepresentationRenderer)
    {
        IRepresentationClass repClass = (geoLayer.Renderer as IRepresentationRenderer).RepresentationClass;
        IRepresentation rep = repClass.GetRepresentation(mapFeature, mapContext);
        if (rep.Shape.GeometryType == esriGeometryType.esriGeometryPolygon)
        {
            IPolygon polygon = (IPolygon)rep.ShapeEdit;
            rep.Shape = feature.Shape; // Geometry for representations just copies from another features
            rep.UpdateFeature();
            rep.Feature.Store();
        }
    }
}


texts: Does anyone know about a plugin that can import a MySQL database into QGIS?


texts: Limiting analysis to particular area of larger DTM using Whitebox GAT?

texts:
I am new at using Whitebox GAT, which I find very interesting.

When one runs an analysis, is there any possibility to limit the processing to a section of the input DTM, or does the latter have to be clipped exactly to the extent of the study area?


texts: Automate converting labels to annotation in ArcMap at multiple scales?

texts: I need to be able to convert labels to Annotation feature classes at multiple scale levels for both address and street data. Similar in concept to generating tiled layers for a web service, I have 5 to 10 pre-defined extents that I need to export annotation for. The resulting annotation will be used in an ArcEngine application that is restricted to those zoom levels. When done manually, I turn on labels for the layers I want labeled, right-click each layer and choose 'Convert Labels to Annotation, save them to an annotation feature class, then repeat for each scale level. Does anyone know of a way of accomplishing this programmatically? All of my code so far has been in Python, but if necessary I'm open to other languages.


texts: Geodetic transformation with Mars data: 'ographic to 'ocentric latitude

texts:
Introduction
I am working with a few raster datasets covering the south polar region of Mars. Some use a coordinate system with aerographic latitudes, which are equivalent in concept to geodetic latitude on Earth (measured relative to the equatorial plane and a line normal to the surface). Others use a coordinate system with aerocentric latitudes, which are like geocentric latitudes (measured relative to the surface and the center of the ellipsoid). In planetary science, these are referred to generically as planetographic and planetocentric coordinate systems.
I would like to apply a geodetic transformation to the 'ographic data so that it will properly align with the 'ocentric data (there can be an offset of several km in the polar regions). However, I have been unable to figure out how to properly define projections for the two systems to use for a transformation.

What I've been able to figure out so far

Typical ways of defining coordinate systems (PROJ, WKT) do not seem to have a way to explicitly say whether latitude is 'ographic or 'ocentric (i.e., geodetic or geocentric).
Defining a geocentric coordinate system (such as with a GEOCCS in WKT) means defining X, Y, and Z coordinates rather than lat/lon as I need to do.
PROJ allows you to convert between geodetic and geocentric latitudes, but it only works on pre-defined ellipsoids for the Earth, not Mars.

Question
How can I run a geodetic transformation of raster data between 'ographic and 'ocentric latitude systems with a custom ellipsoid?
The way I envision this going is:

define a polar stereo projection in which latitude is 'ographic, then assign that to the raster
define a polar stereo projection in which latitude is 'ocentric
use gdalwarp or some other utility to re-project the 'ographic raster to be 'ocentric.

... but maybe something else is needed, too.


texts: qgis2leaf plugin & data defined properties style

texts: I exported a map with data defined properties style from qgis2leaf, but it doesn't appear like the following map. Only dots appears on the openlayers map. Is it possible to visualize this kind of syle on the exported map?


texts: I was trying to create watershed from a DEM that used polylines. So I converted that to raster and the raster to TIN. I delineated the TIN to fill the gaps and then returned it to raster. Then I got the flow direction, flow accumulation and used the snap pour point. I used the flow direction and the snap pour point data to create a watershed and got the pictured watershed. I am relatively new to ArcGIS and can't understand what I'm doing wrong.


texts: Can I add a user with limited privileges to my organization in AGOL

texts:
Is it possible to add a user to an ArcGIS Online organization without giving them access to everything that's shared within that organization?

I'm thinking about contractors from other (small) organizations that do not have their own ArcGIS account, and only need temporary access to a limited set of data.

UPDATE: That user should have read or write access to certain selected elements only, let's say one WebMap. Everything else should be hidden, including items that are shared company wide.


texts: Make Layer from Group Layer Result in ArcGIS ModelBuilder?

texts:
I want to use the ModelBuilder with the following scenario. I use the tool "Trace Geometric Network" in the model, and therefore receive a group layer as the output. I want to use the trace information of ONE layer from this output later in the model. So I have to make a new layer with this ONE layer, but I can't connect the group layer to the "Make Feature Layer" tool.

Would you know how to make a new layer, which is part of a group layer in the ModelBuilder?


texts: Dissolving features to only have multiparts when they are within 100 feet of one another

texts:
I  want to have a selection of parcels that is dissolved by land owner and within 3 miles of a polyline. If the same owner owns parcels on opposite sides of a road (within 100 feet or so), I want the parcel to be a multipart. However, I don't want to have multiparts for parcels with the same owner that are >1 mile of each other.

This is part of a large model that I'm putting together. I'm sure python is the answer but I'm still very green with it. Is there a series of tools that will accomplish this?

I'm using ArcGIS Desktop 10.3 with the Advanced License.


texts: How to assign a value for Slope of Different percent?

texts:
I have slope percent varying form 0 to 25 percent.  Now I want to give a value for each slope percent. For example, if my slope percent is 2 % , I want to assign a value of 10.  As a matter of fact, if it is a vector shapefile, then I can easily add a column for assigning those values in the attribute table. But since it is a raster image, I don't know how to assign those values for every slope percent.
Please tell me how to do it ??!


texts: AddJoin_management does not appear to work?

texts:
I've taken the code snippet below from a larger block of code in order test it and made a script tool with it in ArcMap:

import arcpy
#Set join parameter variables
inFeatures = arcpy.GetParameterAsText(0)
inField = arcpy.GetParameterAsText(1)
joinTable = arcpy.GetParameterAsText(2)
joinField = arcpy.GetParameterAsText(3)
arcpy.MakeFeatureLayer_management(inFeatures,"inFeatures_lyr") #make feature layer from feature class
#Join road feature class to buffer distance table
arcpy.AddJoin_management("inFeatures_lyr", inField, joinTable, joinField, "KEEP_ALL")

In setting up the script tool, I've defined the variable inFeatures as a feature class, inField as a field, joinTable as a table, and joinField as a field.  When I run the tool no errors are generated but no join is created.

Is there a problem with this code snippet?

Is the arcpy.MakeFeatureLayer_management unnecessary or causing problems?


texts: Problem when import vector data from GRASS GIS within RStudio

texts:
I am trying to use GRASS GIS 7.6 within RStudio and import a polygon which is in WGS 84 (EPSG 4326). I have created successfully the location and the mapset in GRASS GIS. When I use the command v.import from RStudio and plot the shapefile I'm getting the image bellow.

While I had to see something like this:

The code I used so far is:

setwd("my_dir")

library(rgrass7)
library(raster)
library(sf)
library(stars)

use_sf()

initGRASS(gisBase = "where GRASS GIS folder is",
          gisDbase = "where my data are",
          location = "spain2",
          mapset = "spain2",
          override = TRUE)

spain = execGRASS("v.import", flags = c ("overwrite", "o"), input = "spain.shp")

plot(spain)

Any thoughts?


texts: How to create template_postgis database?

texts:
I wrote the following code to psql (Sql Shell)

createdb -T template_postgis nyc;

but does not work because syntax error .

trying to make an example of the links related aim GeoServer.
I have postgis 2.0 and postgresql 9.4

http://docs.geoserver.org/stable/en/user/gettingstarted/postgis-quickstart/index.html


texts: Changing "Default" label in label manager of ArcGIS Desktop?

texts:
Does anyone know if there is a way to change the Default Label in the Label Manager for ArcGIS 10?

Otherwise is there a way to mass label features. I have about 500 feature layers in this MXD and in each of the layer tables is a column named "Label".  Because the locations were geocoded, the default uses the text string from the column "Loc_name".  I'd like it to use my label column for all of them.

In the label manager, I have the option to label everyone with the default, thats why i originally asked if that default can be changed.

I've tryed to copy the text string from my column into "Loc_name" but the feild length is set at 14 and it is un-editable and 14 is not enough characters.


texts: Snap points to points of different layer (based on attributes, not nearest)

texts:
I have 2 layer of points, PointLayer1 and PointLayer2. PointLayer1 has to be snapped to PointLayer2.

I have a table that has IDs of PointLayer1 that has to be snapped to PointLayer2.

Sample Table:

ID1        ID2
20         2445
22         2557
23         2456

At the end, all the nodes form PointLayer1 with ID1 has to be snapped to PointLayer2 with ID2.

Completed Code so far:

import pandas
import arcpy

all_snapped_nodes = 'all_snapped_nodes.shp'
all_new_nodes = 'AllnewNodes.shp'

new_old_node_conversion = pandas.read_csv('new_old_node_conversion.csv',index_col = 0, header = None)
#more code thats snaps IDs from layer1 to layer2.


texts: How to reclass a raster with reclassify grid values in QGIS?

texts:
I'm new in QGIS and I'm trying to reclass a raster with reclassify grid values.
I want to reclassify by table, however fixed table has only 3 rows and I need some more, but button to add is not activated.

Do you have any idea about how add more rows?


texts: Change management of aerial photos of same fields

texts: I am looking for an automated solution with capable to stitch aerial photographs and the result can be integrated into Google Maps. Furthermore it has an ability to show discrepancies on a new layer from photos taken  different points in time. For this demands which GIS software could be a solution? Possible to automate this workflow?


texts:
I am trying to perform a slope analysis using Python and GDAL.
The DEMs have a CRS of EPSG:4269 which is NAD83 and a geographic coordinate system. Previous attempts to perform a slope analysis resulted in unexpected results. I am trying to calculate slope as a percent but the results from the slope raster were in the thousands. After some research, I discovered that I should reproject the rasters to a projected coordinate system which I have done in the code below. However, this still produces unexpected results. Some values look great, near zero where I could tell the surface was flat but then other areas would have a slope of 100 or more and that does not seem right to me. I have provided the code I am using to reproject the rasters and perform the slope analysis on.
# reproject raster to projected coordinate system
gdal.Warp(reproject_raster, input_raster, dstSRS='EPSG:32119')

# perform slope analysis
opts = gdal.DEMProcessingOptions(alg='ZevenbergenThorne', slopeFormat='percent')
gdal.DEMProcessing(slope_raster, reproject_raster, 'slope', options=opts)


texts: Unable to retrieve NDVI coloured images to local storage using Google Earth Engine API

texts:
I am currently developing an application which requires  LANDSAT NDVI coloured images to be retrieved  to local storage. to carry out the task i wrote a script in python using earth engine API (source code attached) in which i am able to retrieve images as attached. please note that the retrieved image does not contain NDVI color information (values).

I have two queries with regard to this.

does the API support such queries which get us the NDVI colored data ?
using map.addlayer function i am able to display the coloured image as per requirements. but here i am facing difficulty in getting the file to local storage.

Program:

import ee
import ee.mapclient
import datetime
import os
import wget
import sys
from subprocess import call
ee.Initialize()
ulx = 77
uly = 10
lrx = 75
lry = 8
startyear=2012
endyear = 2013
month = 4
for x in range(2012, 2014):
    startyear=x
    endyear=x+1

    regionname = 'alappuzha'
    bands= ['B2', 'B3', 'B4', 'B5']
    output=regionname+'_'+str(endyear)+'_'+str(month)
    os.chdir('/home')
    region=[ulx,lry], [ulx, uly], [lrx, uly], [lrx, lry]
    strregion=str(list(region))
    def getle8(img):
        return(img.normalizedDifference(['B5', 'B4']))
        #return(img.expression('float(b("first_b50") - b("first_b40")) / (b("first_b50") + b("first_b40"))'))
    unzippedfilename=output+".le8.tif"
    if(os.path.exists(unzippedfilename)):
        sys.exit("File exists:"+output)
    le7 = ee.ImageCollection("LC8_L1T_TOA").filterDate(datetime.datetime(startyear, 11, 1),
                          datetime.datetime(endyear, 11, 1)).map(getle8);
    le7a=le7.mean().select([0], ['le1']).multiply(ee.Image(1000)).int16();

    data=le7a
    print(output+' Processing....      Coords:'+strregion)
    path =le7a.getDownloadUrl({
        'name': output,
        'bands': B2,B3,B4,B5
        'region': strregion
        });
    print("Downloading "+output)
    wget.download(path)
    zipstatus=call("unzip "+output+".zip",shell=True)
    if(zipstatus==9):
        sys.exit("File exists:"+output)
    os.remove(output+".zip")
    print(output+' Finished!')


texts: Field Calculator Python Expression Issues

texts:
I'm pretty new to Python and this is my first post on this forum so I could use some help. I am currently attempting to run some code in the Field Calculator in a model in order to replace string characters (hyphen, period, space) with an underscore ('_'). Here's what I have now.

def customReplace(Layer):
repList = ['-', ' ', '.'];
rep = '_'
for repList in Layer:
    return Layer.replace(repList, rep)
else:
    return Layer

Im getting a syntax error: ERROR 000539: SyntaxError: unexpected EOF while parsing (, line 1)
Failed to execute (Calculate Field (2)).

Anybody out there give me a hand?

By the way, I defined the above function in the pre-logic script code area and I'm calling it via

customReplace( !Layer! )


texts: how to label only number using vbscript in arcgis 10.1

texts:
I have a GPS station feature dataset and I need to label using field NAME. However, I need to label only the number, for example in the field NAME the name are GPS123, STATION 239, STN 325 but the label need only appear 123,239 and 325. I already tried using this VB script expression but its not work properly.

Function FindLabel ([NAME])
  FindLabel = [NAME]
  FindLabel=mid(FindLabel,4)
End Function


texts: Apply Symbology from Layer Not Working in Script

texts:
I am working on applying symbology to multiple layers from lyr files for multiple MXDs within the IDE. I have been able to access the MXDs and iterate through and change the definition of "mxd". I am receiving the following error, however, when the ApplySymbologyFromLayer_management process is run:

ERROR 000732: Input Layer: Dataset FCname does not exist or is not supported

My script is as follows:

# the MXDs to change, accessed and processed individually as 'i' is script not shown
mxd = arcpy.mapping.MapDocument(r'C:\Path\to\mxd\doc.mxd')
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.name == 'FCname':
        arcpy.ApplySymbologyFromLayer_management(
            in_layer="FCname",
            in_symbology_layer="C:\Path\to\lyrs\FCname.lyr")

I have repeated the same ApplySymbologyFromLayer_management tool within ArcMap successfully. Any ideas on why this process is not succeeding within my script?


texts: Does MapInfo have thematic base maps like the templates you get with ArcGIS data and maps?


texts: Python query to find records with more than one alpha character in a field

texts:
Using ArcGIS Desktop 10.2.  I have a field in a table of several thousand records.  The field is made up of alphanumeric characters.  I'd like to perform two queries using python on this field:

Identify the records in this field that contain more than one alpha character
Identify the records in this field that end with more than one alpha character


texts: Changing stacked labels new line character to other than comma using ArcMap

texts:
Trying to execute some stacked labels with little bit of Python from the examples given here:
https://desktop.arcgis.com/en/arcmap/10.3/map/working-with-text/about-building-label-expressions.htm
The example for Python stacked labels uses a comma to denote the new line character but some of my fields have commas in them throwing it off.
Can the character be other than a comma?
Tried just using a blank space and it didn't accept that as an argument. I know semi colon is generally bad Python practice and tried it anyway (didn't work, wouldn't parse).
Suggestions?


texts: Datasets with same geographic and projected coordinate systems, appearing in different locations

texts:
I have two different sets of data (multiple shapefiles in each), when I import dataset 1 lets call it, into ArcMap it displays this error message:

I then use Define Projection tool on dataset 1 which adds its geographic as well as projected coordinate system which is exactly the same as dataset 2 yet they plot in different locations within ArcMap, but they should overlay with each other:

I have also tried to Define a coordinate system with the Define tool followed by the Project tool on dataset 1 but that doesn't work either. Both datasets should be plotting in the same location. Any solution to this?


texts: Troubles replacing null values with specific numbers

texts: After i've used "zonal statistics as table" on some buffers, some values got lost (which seems to be a common problem when buffers overlap as mine does). Now, i would like to replace these Null values by using the "Field calculator". I would like the cells with null values to take on the last value that was not a null value - eg. in the picture attached to this question, i would like the first set of Null values to get the value that preceeded them (in this case 0,655745625495911). I am aware i can do it manually but as the table consists of 330000 rows i would rather not. So, can anyone help me with a pythoncode to change these null values to the value that preceeds them?


texts: Certain fields not show in .csv when using Display XY Data

texts:
I have a .csv file with multiple well locations in Australia. The data is in GDA 1994, although I have put in the easting and northings into GIS manually in UTM Zone 52S and it plots them right where they belong. Among the data columns are three different type of coordinates: lat/long in d m s; easting and northing in meters; and lat/long in decimal degrees. I can manually enter a well's coordinates with the easting and northing and they plot where they should just fine. I have around 3000 wells so manually entering them all will take forever.

Unfortunately, when I bring in the .csv, right click, and go to Display XY Data the easting and northing columns/fields are not available. Both lat-long fiels show up, but of course they plot incorrectly. I have tried different coordinate systems in the base map, doesn't work. I have tried just bringing them in, then projecting them. No surprise, that didn't work. I even tried deleting all the other columns except the easting, northing, and name. Then nothing was available.

Any ideas as to why it would not want the easting and northing?

The map coordinate system is set to UTM, so I don't understand.


texts: HTTP 500 error when uploading shapefile to GeoServer with C#

texts:
I receive an HTTP 500 internal server error when attempting to PUT a shapefile using REST API. I can manually import the file, publish and preview it. So the file is fine.

GeoServer log:

2016-05-24 12:20:25,921 INFO [catalog.rest] - PUT file, mimetype: application/zip
2016-05-24 12:20:27,140 INFO [catalog.rest] - Auto-configuring datastore: testshapes
2016-05-24 12:20:27,140 ERROR [geoserver.rest] -
java.lang.NullPointerException
at org.geoserver.catalog.rest.DataStoreFileResource.updateParameters(DataStoreFileResource.java:515)
at org.geoserver.catalog.rest.DataStoreFileResource.handlePut(DataStoreFileResource.java:271)
at org.restlet.Finder.handle(Finder.java:298)
at org.geoserver.rest.BeanDelegatingRestlet.handle(BeanDelegatingRestlet.java:38)
at org.restlet.Filter.doHandle(Filter.java:105)
at org.restlet.Filter.handle(Filter.java:134)
at org.restlet.Router.handle(Router.java:444)

Fiddler RAW view of request and response:

PUT http://localhost:8080/geoserver/rest/workspaces/test/datastores/testshapes/file.shp HTTP/1.1
Content-Type: application/zip
Authorization: Basic YWRtaW46Z2Vvc2VydmVy
Host: localhost:8080
Content-Length: 205898
Expect: 100-continue
Connection: Keep-Alive

HTTP/1.1 500 Internal Server Error
Transfer-Encoding: chunked
Server: Jetty(6.1.8)

0

C# code:

public bool UploadShapeFile(string workspace, string newStoreName, byte[] fileContents)
{
    byte[] localShapeFile = fileContents;
    String sUrl = GEOSERVER_HOST + "workspaces/" + workspace + "/datastores/" + newStoreName + "/file.shp";
    WebRequest request = WebRequest.Create(sUrl);
    request.ContentType = "application/zip";
    request.Method = "PUT";
    request.Credentials = new NetworkCredential(GEOSERVER_USER, GEOSERVER_PASSWD);

    request.Headers.Add("Authorization", "Basic " +
    Convert.ToBase64String(Encoding.ASCII.GetBytes(GEOSERVER_USER + ":" + GEOSERVER_PASSWD)));

    Stream requestStream = request.GetRequestStream();
    requestStream.Write(localShapeFile, 0, localShapeFile.Length);
    requestStream.Close();

    WebResponse response = request.GetResponse();
    return false;
}

Usage:

var c = new GeoServerClient.Endpoint();
var f =
System.IO.File.ReadAllBytes(@"C:\path\to\file\1-2-2015-34-municipal-shp.zip");
c.UploadShapeFile("test", "testshapes", f);


texts: Referencing layer in Table of Contents using ArcPy?

texts:
I'm looking for the best way to reference a layer in the TOC using arcpy and, if it is possible, that is different than the following:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
layers = arcpy.mapping.ListLayers(df)

Reference layers:

mylayer1 = layers[0]
mylayer2 = layers[1]
...

It comes out that in some cases I may not know the index in the list [0], 1, etc... so I need to reference it using other identifiers.

In previous versions I used to reference the layer first as a file in disk, then add the layer in the TOC and, at last, reference it using the same variable. However, in versions such as 10.2, it does not work and I need to reference in other ways, but the way explained in the example above does not convince me.
Other examples such as the ones explained in the following post are usefull in some cases
(click here for an example but I think there should be a simplest way that does not require loop structures...

So, is there a better way?


texts: Why does GRASS eliminate polygons when using v.generalize?

texts:
I am facing a dilemma while trying to use v.generalize to simplify polygons which are similar in size and structure to US Counties or States.  Often, I lose polygons from the data set after running either a simplification or smoothing method on them.  What troubles me is that there does not seem to be any consistent pattern to this behavior.  In some cases the same set of commands applied to different data sets causes causes completely different results.  As a result, I am unable to develop a repeatable workflow and am reluctantly starting to look for other simplification tools that preserve topololgy.  Here is an example:

# LOAD AND CLEAN
v.in.ogr dsn=counties.shp output=county
v.build.polylines input=county output=county_polylines
v.clean input=county_polylines output=county_clean tool=bpol,rmdupl type=boundary

# SMOOTH AND GENERALIZE
v.generalize input=county_clean output=county_smooth method=sliding_averaging look_ahead=7 slide=1
v.generalize input=county_smooth output=county_general method=douglas threshold=.005

Each of the v.generalize steps results in some of the polygons disappearing.  I'd like to understand why this happens.  Is it due to features being loaded into insufficient memory?  Are my parameters "tight" enough to cause features to be excluded because they've met some sort of geometry criteria? (Note that I'm losing some very large features.)  Is there a problem with the topology of the dropped features which causes them to get erased?  Related to that last question, I assume that build.polylines and clean operations that I do initially should take care of this, but I'm not entirely certain, and I can't seem to find good instructions that address this.  Are there any?

Appreciate any insights people can give, especially those based on personal experience with the v.generalize tool.

Roger


texts: Performance testing ArcGIS Desktop in Citrix environment?

texts:
I have a project (ArcGIS 10.0) for a client that could potentially be deployed in a Citrix environment. It will contain a fairly complex Add-In that will, among other things, import data, allow them to perform reclassifications on rasters, digitize features, and ultimately let them setup, through the Add-In interface, parameters for a 3D modeling scenario in a package outside of ArcMap. These setup parameters will probably be written to a XML file via .NET in the Add-In. The application of course isn't developed yet (as we are in the early stages of the project), but we need to test as much of the future functionality of the application as we can in their Citrix environment and on their local PC installs as well. We have test data to use.

How do I go about performing tests that can provide metrics to compare performance of ArcGIS in their Citrix environment versus their local ArcGIS installs - are the there any software packages and/or methods/protocols I can employ?

We have thought about this some more and are thinking of doing some qualitative tests such as having bookmarks in the ArcMap project and having the user across the country time the draw speed, see how editing of point and line features performs, and 3D display in ArcScene. These are obviously very qualitative and subjective, but it's one approach.


texts:
I am trying to follow this course "https://www.youtube.com/watch?v=IfZSyBf6o84", starting from 20:00, and failed to reproduce the clump function based file, as the warning message shows:
Anyone have idea where the problem is?

The file needed for being processed can be found here:  https://drive.google.com/file/d/1Hu7cAk2cf-EyUQdpZNX8gUmH1JnOqoBi/view?usp=sharing


texts:
I'm trying to use the getMapData function inside an R shiny app, but the function returns NULL whenever I call it.

library(leaflet)

# works
m <- leaflet() %>% addTiles() %>% setView(-80,30,4)
m

# returns null
getMapData(m)

The use in Shiny is no different. To provide an example, this also doesn't work (I've included a couple other outputs to show that the input$map_click is working):

library(shiny)

ui <- fluidPage(
  leafletOutput("map"),
  verbatimTextOutput("mapdata0"),
  verbatimTextOutput("mapdata1"),
  verbatimTextOutput("mapdata2")
)

server <- function(input, output) {

  counter <- reactiveValues(n=0)

  # map as output object
  output$map <- renderLeaflet({
    leaflet() %>%
      addTiles() %>%
      setView(-80,35,4)
  })

  # map as leafletProxy
  m <- leafletProxy("map")

  # store map data on click -- click works, but getMapData not working
  observeEvent(input$map_click,{
    counter$n <- counter$n + 1
    output$mapdata0 <- renderPrint({print(paste0('click worked (click ',counter$n,")"))})
    output$mapdata1 <- renderPrint({getMapData(m)})
    output$mapdata2 <- renderPrint({getMapData(leafletProxy("map"))})
  })

}

shinyApp(ui, server)

Can someone help me figure out what is wrong with the getMapData function in R?


texts: Should the ArcGIS Server Account be in the sysadmin SQL Server Role?

texts:
ArcGIS Server 10.2.2
SQL Server 2008 R2
DBO owned schema

I noticed that our ArcGIS Server account is in the sysadmin server role in our SQL Server 2008.  Does ArcGIS Server place the server account in that role by default?  Seems like too high level of permissions for the server account....

I think a problem I am having is stemming from this.  When I block connections to the database from ArcCatalog, the ArcGIS Server account is still able to connect while I'm doing database management tasks and it causes them to fail.

The screenshot below shows that the database is not accepting connections but the arcgis server account (GIS client machine) is still able to connect.


texts:
if i place the below URL in my browser http://dev.virtualearth.net/REST/v1/Locations?q=bangalore+karnataka+india&key=

The response is below

    {
authenticationResultCode: "ValidCredentials",
brandLogoUri: "http://dev.virtualearth.net/Branding/logo_powered_by.png",
copyright: "Copyright Â© 2014 Microsoft and its suppliers. All rights reserved. This API cannot be accessed and the content and any results may not be used, reproduced or transmitted in any manner without express written permission from Microsoft Corporation.",
resourceSets: [
{
estimatedTotal: 0,
resources: [ ]
}
],
statusCode: 200,
statusDescription: "OK",
traceId: "f54d507047604f1585db1e00254b9a7d|SINM000007|02.00.183.2300|SINMSNVM000034"
}

Not getting the resources

Any issue in the URL?


texts: Is there demand for a service to publish web maps with lots of shapes?

texts:
Is anyone in the GIS community trying to publish choropleth web maps with a large (10s to 100s of thousands) number of polygons? and if so what's your current approach and what are the shortcomings of that approach?

My company is thinking about reviving an old project, but would like to get a sense of the demand before we invest too much time in something that may not be needed.

Features of the service would include:

Support for millions of polygons.
Ultra fast updates (recolor the maps in < 1 second, our goal is 60ms).
Realtime query support (mouseover or click a polygon to display it's info).
Bring your own data (upload your own shapes, classifications, etc)
Built in Census Data (all census layers and data preloaded)
Web access (no coding required)
API access (REST API with clients in Python and Javascript to start)
Overlay (display your published maps as overlays in Google Maps, Bing Maps, Open Layers, etc)

The core of this system is already working, but it would be a pretty major push wrap it in something useable.  We're unsure at this point if there is enough demand to justify the effort. My apologies in advance if this isn't the appropriate forum for this sort of question.


texts: Reading various (non-.shp) vector types using ogrinfo

texts:
I am using the OSGeo4W shell to call ogrinfo on a directory (geodata) to see what vector files it includes:

As you see, it lists two shapefiles and it tells me that it is using the shapefile driver; I assume that is the default driver.

However, my folder includes a KML and a GeoJSON as well.

How would I go about including vector files other than shapefiles in the output as well?


texts: Open iFrame of Links in Map Feature Data in MapBox

texts:
I have a MapBox map of points, with a column called "iframe_link", where there's a different link for each point. I made an iframe that works for one page, but I want to open an iframe for each point.

I'm using the following line of code to try to open the links as an iframe, and it's not working:

'<iframe src="e.features[0].properties.iframe_link"</iframe>'

Any thoughts?


texts: Using Navbar Dropdown Button Checkbox in Mapbox?

texts:
I am creating a map application in the Bootstrap framework. I have a drop-down with 3 buttons to toggle layers on and off. I am using a glyphicon as the checkbox with text to the side. If you click the checkbox, the layer switches as expected. However, if you click on the text it changes the checkbox, but does not control the layer toggle. I have tried wrapping the label with a span but this did not work.

Any suggestions?

HTML:

<li class="dropdown"
BASEMAP

<input type="checkbox" class="glyphicon glyphicon-unchecked" id="mapLabel" value="satellite">SATELLITE

<input type="checkbox" class="glyphicon glyphicon-unchecked" id="mapLabel" value="dark">DARK

JavaScript:

$("#mapLabel").click(function(event) {
console.log("Click");
    layerClicked = window[event.target.value];

        if (map.hasLayer(layerClicked)) {
            map.removeLayer(layerClicked);
        }
        else{
            map.addLayer(layerClicked);
        } ;
});


texts: How to access Google Earth Historical Tiled Imagery Url to other Map APIs or Cesium.js?

texts:
Does anyone know how to access Google Earth Historical Tiled Imagery url?

Then, the users can load the xyz tiled imagery by using different map APIs such as OpenLayers 3, Leaflet, and Google Map APIs or Cesium.js?

I assume that I can make a request of the tiled imagery by the url similar as following: http://mt{1-3}.google.com/vt/lyrs=p&x={x}&y={y}&z={z}&time=yyyy-mm-dd hh:mm

or

http://mt{1-3}.google.com/vt/lyrs=p&x={x}&y={y}&z={z}&begin=yyyy-mm-dd hh:mm&end=yyyy-mm-dd hh:mm

Did anyone try to access Google Earth historical imagery tiles? or it's no way to access it by sending url requests?


texts: Copying and merging lines in one go using ArcGIS for Desktop?

texts:
I want to copy and merge selected features (lines) at a time in destination layer (line).

Is there anyway to do this in ArcGIS 10.2 for Desktop?

Update: My current method involves two clicks using copy and merge.
I used to copied selected linear features and paste to destination layer(selection remains same) and merge it via Editor toolbox.


texts: Aggregating data in 100m*100m blocks on map

texts:
I get customer lat long from an app and my use case is I need to display number of users present in 50m*50m or 100m*100m(deployment time decision) blocks on map.

Some what like below image. No zoom in/zoom out required

So far I have considered MGRS Id (https://en.wikipedia.org/wiki/Military_grid_reference_system)

But It don't provide me 50m*50m block aggregations.
Can someone please suggest me some algo by which I can do this aggregations.


texts: Calculating percentage of points closest to each census tract

texts: Hi if anyone can offer advice on this it would be amazing: I have a layer of geocoded addresses around a central point. I am trying to classify census blocks in the area by the percentage of these addresses that are closer to the centroid of each census block than they are to the central point. Doing it for one census block is easy - is there a way to automate the process for 1000+ blocks?


texts: Mark line which is not connected when intersecting each other

texts:
I wanna create a sign or symbol or anything to inform the map reader that the intersecting line is not connected to other line, eventhough they have same colour.

Here is an image of my work.

I am using QGIS 3.12


texts: arcpy.da.UpdateCursor where clause persists

texts:
I have some unexpected behaviour in an UpdateCursor. I'm iterating over quite a large dataset, and I've found that Python freaks out less if I break up the iteration into chunks:

env.workspace = r"ROADS.gdb"
arcpy.MakeFeatureLayer_management("JUNCTIONS", "jn_layer")

for i in range(-70000, 1500000, 10000):
    strFilter = "UFI >= " + str(i) + " AND UFI < " + str(i+10000)
    print strFilter

    with arcpy.da.UpdateCursor("jn_layer", ["OID@","UFI"], strFilter) as upCursor:
        for rowInt in upCursor:
            # do stuff

I would expect this to create a new cursor for each iteration with different filters, as per the print line:

UFI >= -70000 AND UFI < -60000
UFI >= -60000 AND UFI < -50000
UFI >= -50000 AND UFI < -40000
...etc

The problem is the first iteration works as expected, but each subsequent iteration of i returns no rows - it spits out each subsequent filter value but does nothing inside the with loop. Once it finishes, if I change the initial value:

for i in range(-60000, 1500000, 10000):

then it processes the first 10000 values (which were ignored before), but again the rest of the iterations of the with loop return no values. It feels like it's keeping the first filter on jn_layer and trying to apply the second filter to that filtered set, which of course returns no rows since there's no overlap between the ranges.

Is this the expected behaviour? I would have thought creating a cursor on a layer would only filter the layer for as long as the cursor exists, and the with keyword should be destroying the cursor before the filter changes...? Any advice would be most welcome.


texts:
I have a script that I've used to import information in a CSV to a shapefile.

What I'd like to do is make it so that this script will run on all CSVs in a specific folder location, creating multiple shapefiles.

This is a bit beyond my current programming skills. Can someone point me in the right direction? The current script that imports from a single CSV is below.

import arcpy
import csv
import os
print "arcpy imported"

"""
# Author: John K. Tran
# Contact: jtran20@masonlive.gmu.edu

Input CSV can look like:

Lat,Lon,First Name,Last Name
12.34,56.78,Joe,Smith
0.98,7.65,Jane,Doe
65.43,43.21,Bob,Sagat

Make sure 'Lat' and 'Lon' fields are the 1st and 2nd column in the CSV
respectively.

"""

incsv = arcpy.GetParameterAsText(0) # Change this to the path of your CSV file.
outfc = arcpy.GetParameterAsText(1) + "\\" + arcpy.GetParameterAsText(2) + ".shp" # Change this to the path of your output FC.

spatialref = arcpy.SpatialReference(4326) # Create the spatial reference object as WGS84. Can modify if desired.

if not arcpy.Exists(outfc): # Create the output feature class if needed. (if this is not an existing path, create it)
    arcpy.CreateFeatureclass_management(os.path.dirname(outfc),
os.path.basename(outfc), "POINT", None, None, None, spatialref)

csv.register_dialect("xls", delimiter=",", lineterminator="\n") # Register the dialect for native CSV syntax in Microsoft Excelself.
f = open(incsv, "r")
reader = csv.reader(f, dialect = "xls")

headers = reader.next() # Read the first line as the header names.

for header in headers[2:]: # Add fields for remaining columns if needed. Default is TEXT field.
    arcpy.AddField_management(outfc, header, "TEXT")

cursor = arcpy.da.InsertCursor(outfc, ['SHAPE@XY'] + headers[2:]) # Create
InsertCursor.
count = 0
for row in reader:
    if count % 1000 == 0:
        print "processing row {0}".format(count)
        Ycoord = row[0] # Make sure 'Lat' is in the 1st column.
        Xcoord = row[1] # Make sure 'Lon' is in the 2nd column.
        newrow = [(float(Xcoord), float(Ycoord))] + row[2:]
        cursor.insertRow(newrow) # Insert point in FC for each row in CSV.
        count += 1

del cursor
f.close()


texts: Adding information about adjacent features to data using ArcGIS Desktop?

texts:
I am working with ArcGIS 10.0. As a simple example let's say I have a map of the United States. I want to each row of the data which states have adjacent boundaries with that state - so that in the row for South Carolina there would be new variables adjacent_1 and adjacent_2 that list North Carolina and Georgia (or anything similar to this)

Any thoughts? I know I can select on "features that touch the boundary of the selected feature" but then I'm not sure how to add the information to the data.


texts: How to use ArcPy cursor to define each feature sequentially as current selection of layer?

texts:
I want to implement a script that runs through each feature of a feature class, zooms to the extent of the feature, sets the scale to 1:5.000 and adjacent takes a picture of the scene. The steps are in general very easy. My idea is to use the following code lines:

import arcpy.mapping
myMXD = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(myMXD)[0]
layer = arcpy.mapping.ListLayers(myMXD, "", df)

# Zoom and Scale
df.zoomToSelectedFeatures()
df.scale = 5000

My Question is: How can I run through the whole list of features so that the script will zoom on each feature, set the scale back to 1:5000 and takes a picture? My first idea was to use the Search Cursor to run through the attribute table, but how can I than declare that the feature that is actually "visited" by the cursour is the selected one?

Thanks for any help and suggestions.

Regards


texts:
I have a set of layers that were created using the ImageMosaic plugin. All of them have the TIME dimension available, that we use to display animations on the front-end.

Our idea is to start caching all this layers using GeoWebCache, but we've found that when we send the requests for the tiles using the TIME=, in the response headers we get geowebcache-miss-reason: no parameter filter exists for TIME

Even though that we used the REST interface to pregenerate the cache specifying all the available TIME values. The following is the XML that we use to achieve this:

<?xml version="1.0" encoding="UTF-8"?>
<seedRequest>
  <name>rasterdata_service:__LAYER_NAME__</name>
  <bounds>
    <!--For now, only the bounding box of Catalonia-->
    <coords>
      <double>-2495667.977678598</double>
      <double>-2223677.196231552</double>
      <double>3291070.6104286816</double>
      <double>959189.3312465074</double>
    </coords>
  </bounds>

  <gridSetId>EPSG:900913</gridSetId>
  <zoomStart>6</zoomStart>
  <zoomStop>11</zoomStop>
  <format>image/png</format>

  <!-- type can be seed, reseed, or truncate -->
  <type>seed</type>

  <threadCount>1</threadCount>

  <parameters>
    <entry>
      <!--Generate the cached tiles for the specific TIME value-->
      <string>TIME</string>
      <string>2019-11-29T07:00:00.000Z</string>
    </entry>
  </parameters>
</seedRequest>

Also, in the GWC web interface, if we try to seed the cache manually, we don't see the option for specifying the TIME dimension (Even though it's mentioned in the documentation as one of the possible filters), only to select the Style:

Also the TIME dimension appears correctly when we make a getCapabilties request:

<Dimension name="time" default="2019-11-27T14:00:00Z" units="ISO8601">2019-11-26T16:00:00.000Z,2019-11-26T17:00:00.000Z,2019-11-26T18:00:00.000Z,2019-11-26T19:00:00.000Z,2019-11-26T20:00:00.000Z,2019-11-26T21:00:00.000Z,2019-11-26T22:00:00.000Z,2019-11-26T23:00:00.000Z,2019-11-27T00:00:00.000Z,2019-11-27T01:00:00.000Z,2019-11-27T02:00:00.000Z,2019-11-27T03:00:00.000Z,2019-11-27T04:00:00.000Z,2019-11-27T05:00:00.000Z,2019-11-27T06:00:00.000Z,2019-11-27T07:00:00.000Z,2019-11-27T08:00:00.000Z,2019-11-27T09:00:00.000Z,2019-11-27T11:00:00.000Z,2019-11-27T12:00:00.000Z,2019-11-27T13:00:00.000Z,2019-11-27T14:00:00.000Z</Dimension>

Is there any option within the GeoWebCache that I'm not configuring correctly? Or does the GWC not support the ImageMosaic?

GeoServer version: 2.16.0 running from the official dockerfile (Linux)


texts: Create multiple paths by joining a set of start points in one layer to a single end destination in another layer?

texts:
I need to create a batch of paths so that I can calculate the best potential placement of a development site by counting the number of intersecting features along that path.

I have a grid of points in one layer and a single point in another. I would like to create a line for each of the features in the start layer to the single point in the other layer.

Most of points to line features only take one layer as an input or either take two layers with a single feature.

QGIS on MacOS El Cap

Is there a possible solution to this?


texts:
How do I find the coordinates of the centre / centroid of a raster?

I have tried Zonal Geometry but this does not give me the coordinates.

I need to find the centroids of two rasters and create points and measure the distance and bearings between them. i.e. distance and direction between two raster centroids.

I am using ArcGIS Desktop 10.2.


texts: Google Earth Engine API: How to select an image from an image collection based on NDVI value?

texts:
I want to select the image with lowest median NDVI to represent bare ground cover from Sentinel-2 image collection for my region of interest (<5% cloud).

I have a boundary file pre-loaded on Asset called boundary.

I started with:

var start = ee.Date('2017-05-01'); //Dates of interest
var finish = ee.Date('2018-10-31');
var sentinel = ee.ImageCollection('COPERNICUS/S2')
.filterBounds(boundary)
.filterDate(start, finish)
.filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 5));
//to get imagery with cloud coverage below 5%

//compute NDVI for each image
var NDVI = sentinel.map(
    function(img) {
        return img.normalizedDifference(['B8','B4'])
              .rename('NDVI');
    });
print(NDVI)

// now compute median for each NDVI image
var Dict = NDVI.map(
  function(img){
   return img.set('NDVI',
    img.reduceRegion({
     reducer: ee.Reducer.median(),
     geometry: boundary,
     scale: 10})
   );
   });
print(Dict);

var best = ee.Image(Dict.sort('NDVI').first());
print(best);

It turns out Image (Error) Collection.first: Cannot compare values '{NDVI=0.254702177047319}' (Type<Dictionary<Object>>) and '{NDVI=0.18543783789000937}' (Type<Dictionary<Object>>).

Ideally, I want it to return the index of that image from the sentinel image collection in response to the lower NDVI median so I can select it automatically.


texts: Generate polygons from a set of intersecting lines

texts:
This is a simple and quite common question which has already been asked for different purposes (see this link and this too, for example), here, however, we are looking for not a software package but algorithms that we could try to implement say in Python.

So, as shown below a set of lines are mapped (they are already trimmed, BTW).
Algorithms/ideas to generate polygons (as red ones shown)?


texts: ArcGIS: segmenting buffer based on polygon corners

texts:
I have hundreds of square polygons which I've created a single continuous buffer around the polygons. What I'd like to try and do is try and break this buffer down into segments based on the corners of each polygon so that each polygon has 1 corresponding buffer  zone per side. Is there any way to do this?
I've tried breaking the polygons down into polylines, and then attempted to extend the lines so I could cut the buffer up with those lines, but the extend line tool won't extend the lines the way I imagine.

I'm running ArcGIS 10.2.2, with the Advanced License.


texts: How to cluster external geojson (using leaflet.ajax.js)?

texts:
Understand the method of marker clustering but unable to find the right solution for external geojson. I have tried various cluster plugins but nothing seems to work. See example map here.

var geojsonLayer = new L.GeoJSON.AJAX("https://raw.githubusercontent.com/marinespatialviz/marinespatialviz/master/litter.geojson", {
         pointToLayer: function(feature, yx) {
           var icon = L.icon({
             iconUrl: feature.properties.gx_media_links,
             iconSize: [40, 40]});
             var marker = L.marker(yx, {icon: icon});
             marker.bindTooltip('<p>' + feature.properties.description + '</p>');
             return marker;
         }
       }).addTo(map);
var marker = L.markerClusterGroup();
marker.addLayer(geojsonLayer);
map.addLayer(marker);


texts:
I'm looking to use data driven pages to automate some of my workflow, and I'm almost there for the most part.

I have a layer (Chainage - 200m markers along a watercourse) that I'm trying to perform a Page Definition Query on, but with no success. The chainage layer and index layer both share a field ('system' - sort of an in-house code for the watercourse). However when I set that as my field to query on, I have no chainage features appearing at all. The problem is, I can see that the system number is the definitely the same in the chainage layer and the index layer as it is displayed on my map page via dynamic text.


texts: Obtaining the geospatial complement of a set of polygons to a bounding box in PostGIS

texts:
I want to get the difference between a bounding box (single rectangle) and a collection of polygons (a bunch of complex geometries) within that bounding box.

That is, say my bounding box encompasses the USA, and my polygons are the urban centers polygons, I want to get the geometry corresponding to the regions within my bounding box that are NOT a urban center.

I can easily get the urban centers within my bounding box:

WITH bb_uc AS (
    SELECT urban_centers.geometry as geometry
    FROM urban_centers
    WHERE ST_INTERSECTS(urban_centers.geometry,
                        ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326))
)
SELECT bb_uc.geometry
FROM bb_uc;

This query runs in around 10 seconds in my case, and returns ~100k rows with polygon geometries.

Now, when I try to get the "complement" of this spatial dataset with relation to my bounding box:

WITH bb_uc AS (
    SELECT urban_centers.geometry as geometry
    FROM urban_centers
    WHERE ST_INTERSECTS(urban_centers.geometry,
                        ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326))
)
SELECT ST_DIFFERENCE(ST_MakeEnvelope(-124.7844079, 24.7433195, -66.9513812, 49.3457868, 4326), ST_UNION(bb_uc.geometry)) as geometry
FROM bb_uc

This query takes forever to run (I stopped manually after 3h). I had not imagined that this operation would be that complex (basically, apply a mask to my bounding box)

As an aside, note that I'm using an "unnecessary" WITH statement because I sometimes add a WHERE clause in there on my results.

What can I do to dramatically increase my performance to get to the "complementary" geometry of my intersection? ST_Difference man page mentions not to use geometryCollection, hence my use of ST_Union, which is the operation taking forever. However, I can't think of a way around this.

postgis_version: 2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1


texts: arcpy & TraceGeometricNetwork return ERROR 000584

texts:
I am developing a script to automate the creation of a road network following varying criteria.

I have created a geometric network and I would like to get the connection between several points with the Trace Geometric Network tool. Using the gui, the processing works perfectly and I get the expected result. But using arcpy, when I paste the line got with "Copy As Python Snippet" in the result window, I get the following error message:

>>>arcpy.TraceGeometricNetwork_management("C:/arcout/gruyere.gdb/roadFeatDat/roadNetwork","roadNetwork","C:/arcout/gruyere.gdb/roadFeatDat/points","FIND_CONNECTED","#","#","#","#","#","NO_TRACE_ENDS","NO_TRACE_INDETERMINATE_FLOW","#","#","AS_IS","#","#","#","AS_IS")
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 5110, in TraceGeometricNetwork
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000584: Implementation of this Tool's Validate is invalid.
ERROR 000584: Implementation of this Tool's Validate is invalid.
[...]16 times
Failed to execute (TraceGeometricNetwork).

The extensions enabled are

arcpy.SetProduct('ArcView')
arcpy.CheckOutExtension('Spatial')
arcpy.CheckOutExtension("3D")

The error 000584 should concern a custom developed tool, which TraceGeometricNetwork is not.

Any Idea ?


texts: I'm using CMV as a template for my app, I am using the left and top pane, however I need the app to load with the panes closed, is there any way of doing so?


texts: OpenLayers3: Find the parameters from a WMTS service

texts:
I'm trying to add a WMTS layer in an OpenLayers3 project. I'm able to add this layer in QGIS. Here is the service address: http://www.ngi.be/cartoweb/1.0.0/WMTSCapabilities.xml.

In OpenLayers3, I tried several combination of parameters. The last in date is this one:

var ignLayer = new ol.layer.Tile({
     source: new ol.source.TileWMS({
     url: 'http://www.ngi.be/cartoweb/1.0.0/WMTSCapabilities.xml',
     params: {
          'Layer': 'CartoWeb-Topo',
          'Tileset': '3857',
          },
     format: 'image/png'
     })
 });

Do you think am I missing a mandatory parameter? Or maybe the source TileWMS is not the one to use?

Here are the results of GET queries:


texts: postgis st_buffer + defining coordinate system

texts:
I have a table with 1000 points where the latlng coords are in GCS_North_American_1983. I transformed/changed/set the SRID/coordinate system to 2252 (state plane feet of central michigan) with the below queries.

SELECT ST_SetSRID(geom,2252) from drugpoints;
SELECT UpdateGeometrySRID('drugpoints','geom',2252);
ALTER TABLE drugpoints
   ALTER COLUMN geom
   TYPE Geometry(point, 2252)
   USING ST_Transform(geom,2252);
Update drugpoints set geom= st_transform(geom,2252);
select ST_AsEWKT(geom) from drugpoints;

I am attempting to buffer all these points using the st_buffer with a 1000 foot buffers. when I perform the below query and bring it into QGIS it gives me one enormous circle around all of the points(with 1000 records)

select *, st_buffer(geom, 1000) into drugs2 from drugpoints;

when I pull it into QGIS it looks like this

when I run the buffer tool on the drugpoints in arcgis with the coordinate system of 2252 it correctly buffers each point with a 1000 foot buffers

*running it postgres

My question: why does the drugpoints shapefile lose the coordinate system units of 2252 (which I clearly set/defined/transformed...) when I use the st_buffer function on it? Am I missing a coordinate system query within the buffer query?


texts: Assigning signatures in Erdas Imagine 2016

texts:
I am having a problem in creating new signatures in Signature Editor in ERDAS 2016. Whenever I click on the icon to add signatures a pallette pops up that says:

"ERROR GETTING AOI FROM VIEWER, USE AOI TOOLS IN VIEWER TO CREATE AOIs FOR SIGNATURE EDITOR"

I can't proceed from here and I am new to Erdas so I can't comprehend too.


texts: Adding Distance (KM and Miles) to a Map using ArcGIS 10.1

texts: Using ArcGIS 10.1, how do you add and show the distance in mile / km markings alongside the marked trails(KM and miles) while using the aforementioned 10.1?


texts: Get longitude and latitude from a shapefile with Leaflet in R

texts:
I've been trying to map with leaflet to add addCircles() but my shapefile doesn't seem to have the latitude lat and longitude lng parameters so I got the centroids for each city as is the code below but it doesn't seem to show me the centroid points so I put in a df these values but I can't automate to get the points for n cities.

The data is here.

library(stringr)
library(leaflet)
library(sf)
library(dplyr)
quito = st_read("C:/Users/crist/Downloads/Administraciones Zonales/Administración_Zonal.shp") %>%
  st_simplify(dTolerance = 1000) %>%
  sf::st_transform('+init=epsg:4326')

sectores = read.csv("C:/Users/crist/Downloads/sector.csv", header = T,sep = ";", dec = ",", row.names = 1)
sectores

full_data = inner_join(quito, sectores, by = 'NOMBRE') %>%
  mutate(label_map = sprintf("%s
Valor: %g
",NOMBRE, TARIFA_PROMEDIO_POR_HAB_DISPONIBLE_...) %>% lapply(htmltools::HTML))

bins = c(0, 10, 50, 100, 150,250)
pal_quito <- colorBin("Oranges", domain = full_data$TARIFA_PROMEDIO_POR_HAB_OCUPADA_...,
                      bins = bins)

#I tried to do this because my shp didn't have the longitude and latitude variables

full_data$centroids <- st_transform(full_data, 29101) %>%
  st_centroid() %>%
  st_geometry()

lngt_q = c(-78.41782, -78.67333, -78.4823, -78.60407, -78.36822, -78.50851, -78.56278, -78.3023)
lat_q = c(-0.08668143, -0.2179538, -0.1585809, 0.09029626, -0.4124271, -0.2112893, -0.311081, -0.2025039)

full_data$lngt_q =lngt_q
full_data$lat_q =lat_q

leaflet(data = full_data) %>%
  addTiles() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = full_data,
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              fillOpacity = 0.2,
              label = full_data$label_map,
              fillColor = pal_quito(full_data$TARIFA_PROMEDIO_POR_HAB_OCUPADA_...),
              stroke = T,
              highlightOptions = highlightOptions(
                weight = 5,
                color = "#666666",
                fillOpacity = 0.7)) %>%
  addCircles(full_data$centroids, lng = lngt_q, lat = lat_q, color = "#045a8d", weight = 1, radius =1500,
             fillOpacity = 0.2)

[updated]
There's a way to get the longitude and latitude of my shapefile automatically.


texts: Python Helmert Transformation parameter and residual handling

texts: I'm looking for some starting points to implement an algorithm to calculate Helmert Transformation parameters (4- and/or 7-parameter). I found already some basic implementation e.g. https://github.com/rstulug/Helmert_transformation. Currently I'm not quite sure how to deal with remaining errors after the least square fit.


texts: How to save complex expressions in QGIS Field Calculator?

texts:
I have been debugging some long case expressions in the field calculator and getting frustrated because I have to remember to copy the expression before I test it otherwise I end up having to reenter the expression again, fix the typos, ...

Is there a more efficient way of doing this?  Things like the SQL pane of the the database plugin store past queries.


texts: Are there different EPSG codes for 6 digit, 7 digit and 8 digit UTM Zone32 coordinate systems (CRS)?

texts:
I get data sent to me in various flavours of EPSG:25832. Unfortunately, everyone thinks that their own EPSG:25832 is the correct one. Some send data with the 32 at the front, others with just a 2, and most with the 32 left off and only the False easting of 500,000m

So, which is the correct EPSG:25832? 6 digit, 7 digit or 8 digit? ArcGIS recognises all, but have also various names for them too. E.g in Projection / National Grids / Germany I see 4 different "ETRS 1989 UTM Zone 32N.prj" Files for 6 digits, 7 digits, and 8 digits and one without a suffix which is also just the 6 Digit CRS. Also to make things worse there is also ETRS 1989 UTM Zone N32" (((!!! N32 instead of 32N !!!)))with the 32 at the front again!

In Projected / UTM / Europe/ there is also "ETRS 1989 UTM Zone 32N.prj" which doesn´t begin with 32.

Can anyone standardize this for me?  Has sometime tampered with my CRS Directory in order to confuse me?


texts: Table not found ConvertFeatureClass how to project?

texts:
I get a Table [Streets] not found when I try to convert a shapefile to a geodatabase.

I'm thinking I'll try setting the spatial reference like in this post,
http://forums.esri.com/Thread.asp?c=159&f=1707&t=299127

But I'm not sure what the steps are to get a GeometryDefEdit from my converter code

//input
IWorkspaceName pInShpWorkspaceName = new WorkspaceNameClass() as IWorkspaceName;
pInShpWorkspaceName.PathName = shapeFileDirectory;
pInShpWorkspaceName.WorkspaceFactoryProgID = "esriDataSourcesFile.ShapefileWorkspaceFactory";
IFeatureClassName pFCName = new FeatureClassNameClass() as IFeatureClassName;
IDatasetName pShpDatasetName = pFCName as IDatasetName;
pShpDatasetName.Name = shapeFileName;//this is your input file
pShpDatasetName.WorkspaceName = pInShpWorkspaceName;
IName sourceName = (IName)pFCName;
IFeatureClass sourceFeatureClass = (IFeatureClass)sourceName.Open();
IFields sourceFields = sourceFeatureClass.Fields;

//output
IWorkspaceName pWorkspaceName = new WorkspaceNameClass() as IWorkspaceName;
pWorkspaceName.WorkspaceFactoryProgID = "esriDataSourcesGDB.FileGDBWorkspaceFactory";
pWorkspaceName.PathName = GDBPath; //output file
IFeatureClassName pOutputFC = new FeatureClassNameClass();
IDatasetName pDataSetName = pOutputFC as IDatasetName;
pDataSetName.WorkspaceName = pWorkspaceName;
pDataSetName.Name = "Streets";
IFeatureDataConverter pShpToFc = new FeatureDataConverterClass();

pShpToFc.ConvertFeatureClass(pFCName, null, null, pOutputFC, null, sourceFields, "", 1000, 0);

Also is there anything else I should try ?

I guess it was the 5th Param into ConvertFeatureClass, but unfortunately it didn't make any difference.

*****Update*****

So I did the following operation:

ConvertFieldName(sourceFields, "CITY", "CITY3"); //change column name from city to city3

void ConvertFieldName(IFields fields, string oldname, string newname)
{
  if (string.IsNullOrEmpty(oldname)) return;
  if (newname == oldname) return;
  List<string> fieldnames = GetFieldNames(fields);

  IFieldsEdit editfields = fields as IFieldsEdit;
  int index = editfields.FindField(oldname);
  if (index != -1)
  {
    IField field = editfields.get_Field(index);
    IFieldEdit addresspd = field as IFieldEdit;
    addresspd.Name_2 = newname;
    addresspd.AliasName_2 = newname;
  }
}

And also I try this: (Change column name, and change it back)

ConvertFieldName(sourceFields, "CITY", "CITY3");
ConvertFieldName(sourceFields, "CITY3", "CITY");

Seems insane, I know, but for some reason, doing these operations (either changing or changing back and forth) allowed the GDB file to be created.  I guess the code was cleaning up a resource or something, which allowed for the GDB creation.


texts: What is my distance when viewing a scaled-topographic map?

texts:
A close friend from my company asked me a question about a topographic scaled map (1:25.000) that what was the altitude to a location in km (as you are viewing Earth from space) when he was viewing the map from 1 meter away?

It is a really interesting question, I have been working with scaled maps in the GIS world for a long time but I have never wondered what was my viewing altitude from a topographic map...

If you know any formulae, please share with your explanation...


texts:
Can EasyLoader work with an Azure SQL database?

I've created a valid ODBC DSN to the Azure SQL instance, but when I try and use the EasyLoader to create the "Map Catalog", I get an error "Failed to add user MAPINFO default_database is not supported"


texts: Make a relationship based on points and fences

texts:
I have two models on geodjango, one for tracks and one for fencing, I am trying to match both of them but I couldnt get it working as I spected. My target is check each point if it is inside a fence.

My models are

class Track(models.Model):
    avl = models.ForeignKey('avls.AVLUnit', on_delete=models.CASCADE)
    timestamp = models.DateTimeField()
    latitude = models.DecimalField(max_digits=12, decimal_places=9)
    longitude = models.DecimalField(max_digits=12, decimal_places=9)
    altitude = models.IntegerField()
    speed = models.IntegerField()
    course = models.IntegerField()
    satellites = models.IntegerField()
    extra = models.TextField()
    created = models.DateTimeField(auto_now_add=True)
    address = models.CharField(max_length=512, default="-")

    class Meta:
        verbose_name_plural = 'Tracks'

    def __str__(self):

        return "%s-%s-%s-%s" % (
            str(self.avl),
            str(self.timestamp),
            str(self.latitude),
            str(self.longitude)
        )

class Geofence(models.Model):
    TYPE_GEOZONES = (
        (1, 'LINE'),
        (2, 'POLYGON'),
        (3, 'CIRCLE'),
    )
    avl_resource = models.ForeignKey(
        'AVLResource',
        on_delete=models.CASCADE,
    )
    geofence_id = models.CharField(max_length=8)
    geofence_name = models.CharField(max_length=256, unique=False)
    geofence_description = models.TextField()
    geofence_type = models.IntegerField(choices=TYPE_GEOZONES, default=1)
    geofence_data = models.GeometryField()
    is_active = models.BooleanField(default=True)
    created = models.DateTimeField(auto_now_add=True)
    updated = models.DateTimeField(auto_now=True)

I have a query for make it possible, however I am not realy sure if this is ok, and I would prefer it on a django way

SELECT a.latitude, a.longitude, b.geofence_name FROM tracks_track a, avl_resources_geofence b WHERE ST_Within(ST_SetSRID(ST_MakePoint(a.longitude, a.latitude), 4326), b.geofence_data) = True and a.timestamp::date = date '2018-02-23' order by a.timestamp DESC limit 10;

   latitude    |   longitude   |           geofence_name
---------------+---------------+-----------------------------------
 -14.882885000 | -71.312260000 | GF00
 -14.864581667 | -71.325765000 | GF01
 -14.963245000 | -71.375826667 | GF02
 -14.864731667 | -71.325960000 | GF03
 -14.886156667 | -71.313855000 | GF04
 -14.864593333 | -71.325746667 | GF05
 -14.926816667 | -71.347586667 | GF06
 -14.882885000 | -71.312260000 | GF07
 -14.933970000 | -71.346448333 | GF08
 -14.944961667 | -71.378663333 | GF09

Finally, my question is how to build a query using geodjango or postgis to buld a relationship between this both models / tables or having and optmized way on solve this issue without loosing performance

Best regards


texts: How to create a SpatiaLite Layer in QGIS 3 with Python?

texts:
I have an existing SQLITE Database and I want to create several additional geometry Layers with Python.

Within QGIS 2.x I used the following library

from pyspatialite import dbapi2

According to this issue Pyspatialite in QGIS 3 crash QGIS 2.99 does not have the pyspatialite library anymore. With the sqlite3 library there is obviously no spatialite-extension included.

So now I am stuck. Is there another library I can use?


texts: I'm trying to use an old software that only imports the old style USGS DEM info.  I've found somewhere to download an entire quadrant of which I need about 1 percent.  I have a current LAS file of the area I want.  Is there any way to convert it back to the old style DEM point file format?


texts: Reprojecting MODIS Swath data to WGS84 using GDAL

texts:
Im new to GIS and now working with MODIS04L2 Aerosol Optical Thickness image. I need to reproject the image to other coordinate system such as WGS84 for later work but now Im stuck with it. I have try two ways:

Write my own code using python with gdal lib: I follow the tutorial which using GCPs point as paramater for gdalwarp command. The problem is I can't get the GCPs point. The dataset.GetGCPs() function return a zero-length array of GCPs.
Using programs like HEG tool, MRT tool but both of it have issues. When I open MODIS file (with extension .hdf) with MRT tool I get the error "Error in Module: ReadHdrFile".

I download MODIS from ladsweb.nascom.nasa.gov/data/search.html with parameter: Satellite/Instrument:Terra, Group: Terra Atmosphere Level 2 Product, Products: MOD04L2_Leve 2 Aerrosol and Collection: MODIS Collection 5.1. I want to extract subdataset "mod04:Image_Optical_Depth_Land_And_Ocean" from hdf MODIS file to geotiff which reprojected to WGS84.

The sample hdf file http://sdrv.ms/JpVQNh
The log from gdalinfo:http://sdrv.ms/JpTQok

Anyone here has experiment with this process?

It turns out that the gdal package on my computer is out of date so I cannot get the GCPs via gdalinfo. I updated the gdal to 1.10.0 and now everything is working fine.