texts: Pass coordinates parameters to geoserver pgrouting service from Leaflet wms tile layer

texts:
I've configurated a pgRouting service in Geoserver, following the steps on the documentation it's easy to accomplish.

Reading the documentation, existing examples, etc.. I've come to undestand how the call must be made from the Javascript code:

        var myLayer = L.tileLayer.wms("http://192.168.0.34:8080/geoserver/pgroutingtest/wms", {
            layers: 'pgroutingtest:pgroutingtest',
            format: 'image/png',
            transparent: true,
            version: '1.0.0',
            tiled:true
        });
        map.addLayer(myLayer);

The call is ok, no errors returned, but also no image added to the map.

My doubt is, how can I pass the origin and destination coordinates(x1,y1,x2,y2) as parameters to the pgrouting service in order to calculate the route?

In the example that the pgRouting documentation give us, this is done as follow:

var viewparams = [
  'x1:' + startCoord[0], 'y1:' + startCoord[1],
  'x2:' + destCoord[0], 'y2:' + destCoord[1]
];
params.viewparams = viewparams.join(';');
result = new ol.layer.Image({
  source: new ol.source.ImageWMS({
    url: 'http://localhost:8082/geoserver/pgrouting/wms',
    params: params
  })
});
map.addLayer(result);

But this example, is made with OpenLayers 3. How can achieve this with Leaflet?


texts: What is the easiest way to install a really old Nightly Build of QGIS Master (in Windows)?

texts:
I would like to test which nightly build between QGIS 1.8 and 2.0 caused a certain change in QGIS behavior. I know that there are very convenient exe files of weekly builds available on http://qgis.org/downloads/weekly/ . Yet, their temporal scope is limited. They do not go back  in time enough for my purposes.

The archive at http://download2.osgeo.org/osgeo4w/release/qgis/qgis-dev/ goes back to the very beginnings, but offers only *.tar.bz2 format. Is it possible to use them with osgeo4W "local directory" installation? I have tried with the current osgeo4W version, but the files I downloaded to the local qgis_dev subdirectory were not recognized (osgeo4w continues to offer only the most recent and the previous versions for install). I have downloaded also an old 1.8 osgeo4w version but it does not have the GUI and I do not know what to write in the command line in order to install 1.9 nightly builds.

Well, I really hope there is something more convenient than to compile the source code "manually" according to the methods described on github.


texts: Confirming coordinate system/shapefile is correct?

texts:
I am using Census Shapefiles from NHGIS. I am trying to convert points from it to easy to read latitude and longitude - but am having some trouble. The sample data I'm using for now is data from Alabama.

import geopandas
import random
from shapely.geometry import *
import pyproj
poly  = geopandas.GeoDataFrame.from_file('NewShapes/nhgis0103_shapefile_tl2010_010_block_2010/AL_block_2010.shp')
print poly.crs
p = pyproj.Proj("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
for index, row in poly.iterrows():
    LY=row['geometry'].bounds[1]
    BY=row['geometry'].bounds[3]
    LX=row['geometry'].bounds[0]
    BX=row['geometry'].bounds[2]
    print LY,BY,LX,BX
    lon, lat = p(BX, BY, inverse=True)
    print lon, lat

The output is:

{u'lon_0': -96, u'datum': u'NAD83', u'y_0': 0, u'no_defs': True, u'proj': u'aea', u'x_0': 0, u'units': u'm', u'lat_2': 45.5, u'lat_1': 29.5, u'lat_0': 37.5}
-515926.764688 -511955.403286 856239.343387 860999.442418
-88.180859033 17.9110274036

The latitude appears to be correct - but the longitude is way off. This is an example but I have tried other state shapefiles and other blocks within the data and getting similarly weird results.

Does anyone have an idea of what is going on here?

The bounds (-515926.764688 -511955.403286) do not appear to match any coordinate system I'm aware of.


texts:
Is anyone here using navteq data for a GIS app?

I'd like to work with navteq data and using mapnik-tilecache as a WMS server. The reason is, I use OSM data as well on my other server. And as the two servers need to be interchangeable I have to stick to one software technology.

Importing the .shp files in a PostGIS DB and adjust the .xml config file for the mapnik render engine seems to be the right strategy. Has anyone another idea, or can confirm this route?


texts:
I was wondering what sources are available, or if anyone has water quality data for the Great Lakes and specifically Lake Ontario?

We are interested in mapping Phosphorus, Organic Matter, Nitrogen and Dissolved Oxygen concentrations in Lake Ontario as they relate to algal blooms and hypoxia.


texts: Using rasterio and matplotlib how can I create a window for an NDVI image to fit into?

texts:
I am using Rasterio, Numpy, and Matplotlib. I created a window with the dimensions of 800x600 and plotted the band 5 for a Landsat 8 image, that worked just fine. I also have the nir image so then I calculated the ndvi of the red and nir bands from the Landsat8 image. What I am trying to do is put my ndvi image inside a 800x600 window. How might this be done? So far or the ndvi image it's giving me default dimensions of 7000x7000, I want to create the or use the same window of Window(2000,2000,800,6000)

Here is my code so far:

import rasterio
import matplotlib.pyplot as plt
import numpy as np

nirband = r"LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF"

redband =r"LC08_L1TP_015033_20170822_20170912_01_T1_B4.TIF"

#rasterio.windows.Window(col_off, row_off, width, height)
window = rasterio.windows.Window(2000,2000,800,600)

with rasterio.open(nirband) as src:
    subset = src.read(1, window=window)

plt.figure(figsize=(6,8.5))
plt.imshow(subset)
plt.title(f'Band 5 Subset')
plt.xlabel('Column #')
plt.ylabel('Row #')

rast = rasterio.open(nirband)
rast2 = rasterio.open(redband)
nir = rast.read(1)
red = rast2.read(1)

red = red.astype(float)
nir = nir.astype(float)
np.seterr(divide='ignore', invalid='ignore')

ndvi = np.empty(rast.shape, dtype=rasterio.float32)
check = np.logical_or ( red > 0, nir > 0 )
ndvi = np.where ( check,  (1.0*(nir - red )) / (1.0*( nir + red )),-2 )

plt.figure(figsize=(6,8.5))
plt.imshow(ndvi)
plt.title(f'NDVI')
plt.xlabel('Column #')
plt.ylabel('Row #')

This creates two plots one which is 800x600, the other ranges from 7000x7000, I think this is the default for the image. I am trying to have the same dimensions for my ndvi image as I have for my nir band image.

Here is my attempt at saving the NDVI as a GeoTIFF file:

with rasterio.open("LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF") as src:
    naip_data_ras = src.read()
    naip_meta = src.profile

with rasterio.open('MyExample.tif', 'w',**naip_meta) as dst:
    dst.write_band(1, naip_ndvi, window=window)


texts: How I can fix TileOutOfRange" locator="TILECOLUMN in GeoServer 2.17.2 WMTS in openlayers v6?

texts:
I am going to use GeoServer WMTS in openlayers and then I found an example in Openlayers website, It seems everything is correct but I get this error:
<ExceptionReport xmlns="http://www.opengis.net/ows/1.1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 http://geowebcache.org/schema/ows/1.1.0/owsExceptionReport.xsd">
<Exception exceptionCode="TileOutOfRange" locator="TILECOLUMN">
<ExceptionText>Column 21 is out of range, min: 39 max:43</ExceptionText>
</Exception>
</ExceptionReport>

Where should I change in my code, My map zone layers is UTM Zone 38N untill UTM Zone 42N, I chanage EPSG based on UTM zone area but I got 400 bad request again.
const projection = ol.proj.get('EPSG:4326');
const projectionExtent = projection.getExtent();
const size = ol.extent.getWidth(projectionExtent) / 256;
const resolutions = new Array(14);
const matrixIds = new Array(14);
for (let z = 0; z < 14; ++z) {
  // generate resolutions and matrixIds arrays for this WMTS
  resolutions[z] = size / Math.pow(2, z);
  matrixIds[z] = "EPSG:4326:" + z;
}

var wmtsLayers=[
    new ol.layer.Tile({
      title:`${layers_}`,
      opacity: 0.7,
      source: new ol.source.WMTS({
        attributions: '...',
        url: `http://localhost:8080/geoserver/gwc/service/wmts?request=getCapabilities`,
        layer: `${geoServerWs}:${layers_}`,
        matrixSet: 'EPSG:4326',
        format: 'image/png',
        serverType:'geoserver',
        TileLoadFunction:function(tile, src) {
            var xhr = new XMLHttpRequest();
            xhr.responseType = 'blob';
            xhr.open('GET', src);
            xhr.setRequestHeader('Authorization', 'Basic ' + window.btoa(`${geoUserName}` + ':' + `${geoPass}`));
            xhr.onload = function() {
                var objectURL = URL.createObjectURL(xhr.response);
                tile.getImage().onload = function() {
                    URL.revokeObjectURL(objectURL);

    };
    tile.getImage().src = objectURL;
  };
  xhr.send();
},
        projection: projection,
        tileGrid: new ol.tilegrid.WMTS({
          origin:ol.extent.getTopLeft(projectionExtent),
          resolutions: resolutions,
          matrixIds: matrixIds
        })
      })
    })
  ];
 var wmsLayers=[
     new ol.layer.Tile({
         title:`${layers_}`,
           visible: true,
      source: new ol.source.TileWMS({
       url: `http://localhost:8080/geoserver/${geoServerWs}/wms`,
       params: {'LAYERS': `${geoServerWs}:${layers_}`},
 ratio: 1,
        serverType: 'geoserver',

        imageLoadFunction: function(image, src) {
            var xhr = new XMLHttpRequest();
            xhr.responseType = 'blob';
            xhr.open('GET', src);
            xhr.setRequestHeader('Authorization', 'Basic ' + window.btoa(`${geoUserName}` + ':' + `${geoPass}`));
            xhr.onload = function() {
                var objectURL = URL.createObjectURL(xhr.response);
                image.getImage().onload = function() {
                    URL.revokeObjectURL(objectURL);

    };
    image.getImage().src = objectURL;
  };
  xhr.send();
},

      })
   })
];

// The Map
    var map = new ol.Map({
        target: 'map',
       view: new ol.View({
             zoom: 5,
            center:ol.proj.fromLonLat([54.233527,32.343721])
         }),

        layers: [baselayers,],
       overlays: [overlay],

    });

*GeoServerWm is name of worksapce.


texts: How to convert multiple Geospatial PDFs to map tiles

texts:
What is the best way to convert Geo-referenced PDF files (like the Canadian BC http://pub.data.gov.bc.ca/datasets/177864/pdf/ or USGS maps) to a collection of tiled XYZ web maps.

I would like to retain as much detail from the PDF's as possible (ie: sharp lines at say map tile zoom of 18) so an ultimate solution would convert directly from the geospatial PDFs to tiles, however I assume this cant be done due to map stitching steps?

I wish to do this with open-source software where possible (I already have gdal installed).

Please also include recommended map formats for any intermediary steps required.


texts: What's a robust strategy for displaying 3D map data for navigation?

texts:
3D map data, specifically the kind you see on GPS units, or more recently on smart phones, almost always features a 'HUD' like view, which I'm trying to implement in an application of mine. I have OpenStreetMaps data (offline) available to my application, where all the geometry is specified in latitude and longitude.

With 2D maps, there's a 'standard', where almost everyone uses N*N tiles of Mercator projected data. But with 3D data,

I don't know what projection I should be using (I still have to map the latitude, longitude and altitude to x,y,z)
I'm not too sure what's the best way to make data continually available in 'chunks' to a display engine as a user pans and zooms the map view

If it matters, I'm using C++ with the OpenSceneGraph library (which is basically a scene graph wrapper/toolkit for OpenGL). I'd appreciate any advice!


texts: Implementaion of GSCoverageEncoder for WorldImage in Java

texts:
I want to upload a raster image to GeoServer as WorldImage from my Java application. I'm using geoserver-manager v1.4.2 - a Java-based implementation of GeoServer's REST API. It provides convenient methods to work with the  GeoServer REST API.

The problem is that the support for the raster image types is not implemented (only Image Mosaic is supported , which doesn't suite me).

Does anyone has an example implementation of GSCoverageEncoder for an image with world file?


texts: Using Geographic Transformation in Projection?

texts:
I copied my coordinate files into Coordinate folder inside Users\myuser\AppData\Roaming\Esri\Desktop 10.2\ArcMap\Coordinate Systems\....

The problem arises when I try to project using a custom made Geographic Transformation. When I select the Output Coordinate then in Geographic Transformation (optional) there comes nothing.

Where are the transformation parameters created?

Amazingly, this only happened in one of my laptop. Otherwise, this was not the problem before.


texts: Printing the mean NPP value from MODIS imagery in a specific region?

texts:
I am getting started with Google Earth Engine, trying to do something simple: calculate the mean Net Primary Productivity in a region that is defined by a polygon. NPP values are taken from the MODIS satellite data. I can get the MODIS data, select the NPP band, filter the area by a polygon, but I am unable to do the simple act of printing out the mean value of NPP in the region.

Instead of printing out single numeric value (the mean), the Console prints out YAML metadata for an Image (with id, crs, crs_transform, data_type and so forth). I don't want the damn image metadata, just a simple, single, computed value. My JavaScript code is below:

//Get the MODIS data
var mod_npp = ee.ImageCollection('MODIS/055/MOD17A3');
var mod_npp_only = mod_npp.select(['Npp']);

//Define the area of interest
var lmnp_terr =
    ee.Geometry.Polygon(
        [[[35.626418548030756, -3.3142503419625107],
          [35.626418548030756, -3.9090722349373133],
          [36.049392180843256, -3.9090722349373133],
          [36.049392180843256, -3.3142503419625107]]], null, false);

//Filter the MODIS data to the area of interest
var lmnp_npp = mod_npp_only.filterBounds(lmnp_terr);

// (try to) calculate the mean NPP in the filtered area
var mean_npp_lmnp = lmnp_npp.reduce(ee.Reducer.mean());

// (try to) print out a single integer value of mean NPP in the polygon. This fails.
print(mean_npp_lmnp);


texts: MODE_SNAPSHOT or MODE_ONDEMAND? Issues with both, which one to focus on?

texts:
Our web app uses the ESRI JS API. With SNAPSHOT, our feature layers take a long time to show up on the map, especially after using setDefinitionExpression to restrict the feature layer to a date range. If we change it to ONDEMAND, our feature layers show up very quickly, but not reliably. It often takes wiggling the mouse (panning a small amount of pixels) to get the feature layers to show up.

Does this sound familiar to anyone? Which of the two modes would be a better starting point to get things working right? SNAPSHOT got much quicker after we turned off TimeInfo in our feature layers, but it's still noticeably slower than ONDEMAND.


texts: Looking for a simple method to find the shortest path with Shapefiles

texts: I am having shapefiles for a particular map location, the shapefile consists of points and polylines. Now I want to find the shortest path between the points. What is perfect and simple method to find shortest path.


texts: Using gdal2xyz.py to get Raster's pixel center coordinate with it's value

texts:
I am using gdal2xyz to get pixel values translated to a point layer. my values have decimals and when I use gdal2xyz the decimals are removed. so value that were 0.345 become 345 and values that were 1.234 become 1234. I am using it like this:
os.system('gdal2xyz.py -csv ' + Folder_2 + '/Water_Vapor.tif ' + Output+ '/merged.csv')


texts: Removing unwanted characters from coordinate field

texts:
I have some data that arrives via an Excel XLSX file that I import into a feature class.  The Excel file has a field GPS where the users have typed in coordinates, eg 32.1357, -116.81828, which I split into two separate fields in order to use the XYTableToPoint tool to generate the feature class.
The problem is that they're typing the coordinates into the field, so sometimes my coordinate has no space,  32.1357,-116.81828, sometimes it has no comma, 32.1357 -116.81828, etc.  Really I guess it could be anything, until I can provide a solution that eliminates this user error potential (not an option at this point unfortunately).
Current within my UpdateCursor I'm using a series of str.replace() to try to eliminate the potential issues, but I believe there must be a more efficient way to do it.
with arcpy.da.UpdateCursor(in_mem,["GPS","GPS_X","GPS_Y"]) as cursor:
    for gps, gps_x, gps_y in cursor:
        m_gps_y, m_gps_x = gps.replace(", ","#").replace(",","#").replace(" ","#").split("#")
        gps_x = m_gps_x
        gps_y = m_gps_y
        cursor.updateRow([gps,gps_x,gps_y])

How can I remove these unwanted extra characters from a field of coordinates, or extract just the coordinate values for my two new fields?


texts:
How do I view the gridded data provided by gSSURGO in QGIS?

Following the instructions here (File Geodatabase (*.gdb) support in QGIS?), I'm able to import the vector data and tables with no problem but can't figure out how to locate or import the gridded (raster) layer.

Can this be done?

I am running QGIS 2.8.6-Wien on Ubuntu LTS 16.04.


texts: Batch conversion of Shapefiles to Topojson

texts:
Can i convert more than one shapefile .shp to topojson. Can this be performed using batch processing?

I can convert shape file to topojson one by one but i have about 6000 files and i don't want to merge them i want them to be separately converted to topojson.

I can also batch process ogr2ogr but don't know how to perform that on topojson. Is there any way i could perform this task using a folder of is there any alternate solution to this problem.


texts: Changing Style of WMS dynamically with a select box in Leaflet?

texts:
I have a Leaflet WMS Layer
mylayer = new L.tileLayer.wms('https://{s}.url.com/wms?', {
    mylayers: 'Workspace:mylayer',
    format: 'image/png',
    transparent: true,
    styles: 'sld1',
});
mylayer.addTo(map);

I have others styles sld2, sld3, sld4 how can change styles dynamically with a select box. I tried to use .setParams but it does not work:
<select id="styles" onchange="layer.setParams({styles: this.value});">
    <option selected="selected"  value="sld1">sld1</option>
    <option value="sld2">sld2</option>
</select>

This is the code:


texts: Merge Multi-Band Geotiff with different datatypes (with gdal)

texts:
I have three geotiffs:

orthomosaic.tif: a multispectral orthomosaic with three channels:

Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
Band 4 Block=512x512 Type=Byte, ColorInterp=Undefined

DTM.tif (Digital Terrain Model) and DSM.tif (Digital Surface Model) of the same extent with one band each:

Band 1 Block=5001x1 Type=UInt16, ColorInterp=Gray

I want to merge those three into one single file (to consequently retile them and work with smaller tiles having all six channels). This works fine with this command where the parameter -seperate preserves the bands and -o stands for "output":

gdal_merge.py -o merged_file.tif -separate  orthomosaic.tif  DTM.tif  DSM.tif

Output is: merged_file.tif

Band 1 Block=10002x1 Type=Byte, ColorInterp=Gray
Band 2 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 3 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 4 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 5 Block=10002x1 Type=Byte, ColorInterp=Undefined
Band 6 Block=10002x1 Type=Byte, ColorInterp=Undefined

Problem is: The script converts all channels to Type=Byte, which makes the DTM/DSM channels Type=UInt16 useless for me.

Question is: Is it possible to have different data types in different bands in a single geotiff file? And how do I do this when merging several geotiffs?


texts: Using QGIS Processing Modeler to split polygons with lines

texts:
I am trying to use the QGIS Processing Modeler to split polygons with lines.
Here is the workflow I have so far.  It is sucessful.

I then add "Grid" and "Split With Lines" to split the polygons.  No matter what combination of options I select it always fails.  The end result I get for the split polygons is grid of lines.  I have also tried to change the "Parent Algorithms" option to a number of combinations with no luck.

Here is the end product I would like to achieve.

Here is what I think should succeed, but does not.

Here is the result.


texts: Forest Coverage - European free GIS datasets

texts:
Where can I find GIS data for the Forest coverage for the entire Europe? Is there such a product available for free?

I am looking for raster or vector (polygons) with the information of the % of coverage by raster cell or polygon.


texts: Info template displayed far from mouse click

texts:
I've consistently noticed in my ESRI webmap that the info template is displayed below the mouse click. For example, if I click on a polygon the one below it is displayed. Sometimes, I cannot get the info template for a feature to display. Here is a screenshot of the issue, the white arrow shows where I clicked, and as displayed, the feature 2 polygons away is highlighted.

Has anyone ever seen this before, or know how to fix?

I'm relatively new to javascript so not sure what to try. What i've done so far is to search for questions similar to mine. My best guess is that maybe something is messing with the mouseY position on click?


texts: Exporting data from Collect Values output in ArcGIS ModelBuilder?

texts:
I'm currently playing with ModelBuilder. I've added a picture of the model I have at the moment.

This model is currently iterating 6 shapefiles, so each of the 'output values' contain 6 figures in a list. I can't find a way to extract the values from these lists into a table/text file or similar.

Is there a way to do this?

This is what is shown when I open 'output values' after running the model:

I just want to grab those 6 numbers somehow...


texts: Setting DATA string on LayerObj in MapFile by Using Mapserver's PHP Mapscript

texts:
I'm manipulating mapfiles through mapscript. Setting layer attribute values like the opacity works like this:

$layer = $map->getLayer('layername');
//set an attribute value
$layer->set('opacity', 80);

But now I want to set the data attribute as well, I assumed it functions like this:

$layer = $map->getLayer('layername');
//set an attribute value
$layer->set('data', '/path/to/my/shape/file.shp');

But nothing happens, i.e. an image is constructed but it's empty because no data was inserted.

How can I set the data string?

I'm using Mapserver 5.0.3 + Mapscript 4.10

I used the map.getLayer() method instead of the getLayerByName method. While testing on the opacity attribute, I had a mapfile with only one layer, so the getLayer('loadlines') did return a (but in this case the last and correct) layer.


texts: Index out of range error when cleaning up variables

texts:
I'm running a script that parses points from a flat database structure into individual feature classes with related tables.

A search cursor runs through and stores each field's data in a variable that I use to create points with using an insert cursor.  After each loop, I'd like to clean up the variables so that data from the last loop can't be used in each loop thereafter.

I've got an array that I've hard coded with each variable name that I loop through to delete each one.  I've checked these var names 100 times to make sure they're correct (even copied them from where they're declared and pasted them into the array).

I'm getting an error:

  IndexError: list index out of range

I'm using len() to find the length of the array, so I can't imagine how it could be out of range.

    # Clean-up Variables
    cleanup = [ptGeom, ptName, ptType, ptDate, active, ptComm, fdDate, fdCheck, surType,
               fdType, special, lure1, lure2, fdComm, photo, signType, signAmt, platSet,
               platSign, trapType, trapAct, trapStat, capSpec, capSex, capType, trapSet]
    length = len(cleanup) - 1

    for i in range(0, length):
        print i
        if cleanup[i]:
            del cleanup[i]

The error is after it prints 19 (trapType) -- which is a null value according to the attribute table.  The loop isn't having trouble with any of the other null values that it encounters before this line.


texts: Creating feature-linked annotation as either single-part or multi-part

texts: I have two feature-linked annotation features in a geodatabase. Both appear to be identical in every way, except that one creates annotations as single-part and the other creates them as multi-part. I know I could convert one to the other, but I'm looking for an option to create a new annotation feature as either single-part or multi-part. If there is no such option, how does ArcMap decide to which to create an annotation as?


texts: PostGIS - Raster - update all values - what is prefered

texts:
I have added the Raster support to my PostGIS database program. Now, I need to update all values in the raster (one band only) several times a day.

I have found and tested these two options

Use ST_SetValues (http://postgis.net/docs/manual-2.1/RT_ST_SetValues.html)
Use WKB and do UPDATE table SET rast = 'wkb_data_string'; (code extracted from http://www.postgis.org/documentation/raster-doxygen/dd/d23/rt__api_8c-source.html)

Both methods works fine with same results, speeds according to my tests seems also almost the same (but I may have done some errors there).

What is a better way from those twos, or is there some third option?


texts: Editing non-nullable fields of ArcSDE layers in ArcMap

texts:
I'm using ArcSDE layer in my map project in ArcMap 10. Some of the fields of this layer allow null and some don't. I can't edit the non-nullable ones, neither programmatically (IField.Editable is false) nor via Attributes panel.

Is this the restriction of SDE layers or is it possible to tweak some settings/make some changes in DB schema to allow editing of fields that can't contain null values?


texts: Replace attribute value with non NULL values from a joined attribute column

texts:
I have a Geodatabase with a column for pipe diameter. I have a secondary data set with field verified data and a column for pipe diameter. I would like to replace the pipe diameter in the Geodatabase with that in the secondary data set using the Python Parser in the Field Calculator. I have joined the secondary data set to the geodatabase and have many NUll values associated with pipes that have not yet been field verified.

Through online research, without finding anything specific, I have hodgepodged the following:

def IgnoreNull (GravitySewer.Diameter, assets.Diameter)
 if assets.Diameter is None:
      return GravitySewer.Diameter
 else:
      return assets.Diameter


texts:
Does anyone of you know a QGIS-Plugin for Versioning PostGIS-Data?
I know there was pgVersion of Horst Düster but as far as I know this was not updated for QGIS 2.x (at least not for the "normal" QGIS aside from QGIS Enterprise.
The searched functionality is:
Setting up versioning for a Layer
Managing conflicts that were detected (resulting from concurrent editing)

old Plugin pgVersion: https://plugins.qgis.org/plugins/pgversion/


texts: How can I calculate the percent of a raster cell covered by a polygon?

texts:
First, I am going to try and intersect a raster and a polygon shapefile. I'm not 100% sure how to do this (or if I need to do this for the next step), so if you have any suggestions please let me know! Then, my next task is as follows:

If the overlayed/intersected polygon covers > 50% of a raster grid cell, then the grid cell gets a value of "1". If not, the grid cell gets a value of "0".

Does anyone know how I can execute this task in ArcMap 10.3.1?


texts: What are the most tangible benefits to getting your GISP or ESRI Certs?

texts:
Possible Duplicate:
  GISP - Value of being certified?

Looks like I'm only a few presentations or other "contributions to the Profession" points away from being able to apply for my GISP. Has anybody found any actual benefits to getting the certification?

Seems like it might be somewhat arbitrary. Of course it lets people know you have a high level understanding of GIS but beyond that are there any benefits? Salary bumps? Other?

ADDED

Additionally, what, if any are the benefits to getting the new ESRI Certifications? ESRI is obviously the most well known GIS software, but there are others... so if one is looking to move forward in their GIS career, what is the best way to do it as far as certs/continuing education go?


texts:
Suppose we plot a number of points (latitude, longitude) on the google map using Plot Google map package of R. These points are classified on the basis of unique vehicle-id i.e. a point to represent where a particular vehicle stopped. So, the data consists of numerous vehicles. Now, I just wanted to know whether its possible if we click on a particular legend(marker) on the Google map, then all the legend(marker) of that particular vehicle-id gets highlighted.
I just wanted to know whether it could be done or how to do it.

vehicle_id  longt   latit
19969     86.3607   23.80903333
19969     86.346225 23.7921
23384   86.0130875  23.78897883
23384   86.04025    23.840275
23384   86.1597922  23.914832
23384   86.72869833 23.8198765
23395   86.00305    23.766375
23395   86.58922857 23.8085
23395   86.2883     23.82095
23400   86.01176667 23.78473333


texts: Summing pixel values within polygon zone in ArcGIS Desktop?

texts:
My question is the same as Summing pixel values within polygon zone in QGIS? but I'm working with ArcGIS and not QGIS.

I have a raster data set with values of type float.

And a polygon shapefile

I want to calculate mean of all raster cell values that are completely within the shapefile

In order to do so, I have used Extract by mask tool and now I have the following raster with NULL values for cells that are outside the polygon layer.

But there are 3 problems:

All parts of each masked pixel is not within the shapefile.

I want to mask just the pixels that are completely within polygon shapefile.

I don't know how to sum up all cell values in a masked raster?
I don't know how to count the number of such cells and calculate the mean?


texts: Get value from point layer based on start/end of line layer

texts:
I need to get a value from a point table based on the end/start point of a line. This is similar to my question QGIS Field Calculator Aggregation with a intersect query

aggregate( layer:= 'Point_Clip', aggregate:='concatenate', expression:=to_string(RL), concatenator:=', ', filter:=intersects($geometry, geometry(@parent)) )

This doesn't always put the start point (based on the direction of the line) and then the end point.

If I get the start/end points using geom_to_wkt(start_point($geometry))

Can I use this to get the value...I was trying something like

dbquery('Point_Clip','Name','start_point(@geometry) intersects $geometry')

But this gives errors I can't work out. In the end I need it to run in the QGIS Attributes form so that it also fills values from the points at the start/end of existing and new lines using something like

if(regexp_substr(aggregate(
layer:= 'Point_Clip',
aggregate:='concatenate',
expression:= to_string(Name),
concatenator:=', ',
filter:=intersects($geometry,
geometry(@parent))),'(.*),')='',geomintersects('Point_Clip','Name'),regexp_substr(aggregate(
layer:= 'Point_Clip',
aggregate:='concatenate',
expression:= to_string(Name),
concatenator:=', ',
filter:=intersects($geometry,
geometry(@parent))),'(.*),'))

The test data is only about 20 lines and points - it can be downloaded from https://drive.google.com/file/d/12rUV_pAaevs9vATbIsgc9QJvvziRygGN/view?usp=sharing


texts: Analyzing date and time in ArcGIS for Desktop?

texts: I am have a Date/time field (MM/DD/YYYY hh:mm:ss format) in ArcMap 10.2.2 that I want to filter through and show the specific hour periods for accidents. The end result needs to help to show the different time periods of the day that accidents are more frequent. Example: AM vs PM or 6:00AM to 6:00PM vs 6:00PM to 6:00AM. Is there a specific tool in the Toolbox I can use? Or does anyone know a process to help me show these differences on a map.


texts:
I was fortunate to find an SVG of this projection at http://bl.ocks.org/mbostock/raw/4476487/

However, there appears to be little reference to this projection outside of wikipedia. Does anybody know what the Proj4 string would be?


texts: In ArcMap 10.4.1 I would like to make a raster with a simple Gaussian distribution away from a central point - as in the image below. It does not need to be 3D, simply higher values at the centre and lower values to the outside. This will act as a null model to test against the actual distribution I found. The 'Create Normal Raster' function in Spatial Analyst appears to make a raster with randomly distributed values that over-all provide a normal distribution, but I want a surface that represents continuous dispersal away from a point.


texts: How to create a printable map of labelled markers with OSM background?

texts:
I want to create an image of an OpenStreetMap map suitable for printing on A4 at 300dpi.

The map needs to have half a dozen markers/points, and each should be labeled with the name of the point.

I know the name, latitude and longitude of each point.

How can I create such a map in QGIS? Ideally it will also include the OSM credits along an edge. QGIS isn't an absolute requirement: any other free/open source tool, website or script is fine; first choice for Mac, second choice Windows.


texts: Adding graticule overlay to satellite image

texts:
I'm developing a python script to turn Sentinel-1 satellite images in .SAFE-format into easy-to-read jpeg images. The calibration and projection works fine using gdal, but I am having trouble with the last part.

I want the generated image to have a graticule overlay, making it possible to reference geographical spots on the image. I feel that such a normal problem would have a built-in function or some library to solve, but I can't seem to find any. Does anyone have any ideas or tips on this matter?

Since the tif-file have the coordinates of each pixel it would be possible to make a transparent copy of the image, drawing the lat/long lines "by hand" using PIL and then placing this image on top of the other, but this seems like an awful lot of work.

I'm fairly new with gdal.

Using mkgraticule.py i've managed to make the grid lines I wanted. However, they won't help that much if they're not numbered (lat/long for each line). From what I can see mkgraticule don't have this option so I suppose I'll have to find another way.


texts: How to use GDAL to write lat,lon,z from NetCDF to GeoTIFF?

texts:
I have the following NetCDF file:
Format:
           netcdf4_classic
Global Attributes:
           history     = 'FERRET V7.43 (optimized)  9-Jun-20FERRET V7.43 (optimized) 19-Jun-20FERRET V7.43 (optimized) 18-Sep-20'
           Conventions = 'CF-1.6'
           Projection  = Stereographic Oblique Projection with 70.5°N, 40°W as centre
Dimensions:
           X12_251   = 240
           Y20_465   = 446
           SECTOR    = 2
           TIME      = 365   (UNLIMITED)
           SECTOR1_1 = 1
           ATMLAY3_3 = 1
Variables:
    X12_251
           Size:       240x1
           Dimensions: X12_251
           Datatype:   single
           Attributes:
                       units         = 'km'
                       long_name     = 'x'
                       point_spacing = 'even'
                       axis          = 'X'
    Y20_465
           Size:       446x1
           Dimensions: Y20_465
           Datatype:   single
           Attributes:
                       units         = 'km'
                       long_name     = 'y'
                       point_spacing = 'even'
                       axis          = 'Y'
    LON
           Size:       240x446
           Dimensions: X12_251,Y20_465
           Datatype:   single
           Attributes:
                       units     = 'degrees'
                       long_name = 'Longitude'
                       history   = 'From ICE.p42.2019.01.01-05'
    LAT
           Size:       240x446
           Dimensions: X12_251,Y20_465
           Datatype:   single
           Attributes:
                       units     = 'degrees'
                       long_name = 'Latitude'
                       history   = 'From ICE.p42.2019.01.01-05'
    SMB
          Size:       240x446x1x365
          Dimensions: X12_251,Y20_465,SECTOR1_1,TIME
          Datatype:   single
          Attributes:
                      missing_value = -9.999999790214768e+33
                      _FillValue    = -9.999999790214768e+33
                      long_name     = 'Surface Mass Balance (SMB~SF+RF-RU-SU-SW)'
                      units         = 'mmWE/day'
                      history       = 'From ICE.p40.2018.01.01-05'

The file has lat,lon coordinates as 2D arrays (240x446) but are projected to a Stereographic Oblique Projection with 70.5°N, 40°W as centre. The SMB variable is shown on this grid but on a daily basis throughout the year (240x446x365). How can I use GDAL to write the yearly average SMB variable of this file into a GeoTIFF with NSDIC Sea Ice North Polar Stereographic Projection (EPSG:3413)?
I share a Google Drive link with an example file:
https://drive.google.com/file/d/1iPICsNsZ_fLGHuO0gklc8p8ypWaMo-ZO/view?usp=sharing


texts: QGIS 3.4 Add vector styles and attribute fields in a model

texts:
I made a model, which works like this: from a polygon, it extract its edges coordinates, turns it into lines and add an extra line.shp who will store the distance between the points.

Because my skills are basic, I can't figure out, firstly, how to add a new field in the attribute table "distance", the command $length and secondly, how to add style in my new layers.

I have .qpt files with styles stored in my PC. I read the QGIS manual, I load them in the model which does normally all the work, but the result are some default style layers.

Generally, the model, does the job that I want, but I have to add the style to the new layers and add the length field... I want to improve it, to work faster.

Do you have any suggestions?


texts: ArcMap not respecting advanced sort in related tables?

texts: ArcMap 10.2.2 doesn’t respect the advanced sort in related tables, I couldn’t figure out why the ArcMap doesn’t respect the advanced sorted values when the sorted table is related to other feature calss?


texts: Vector to Raster conversion doesn't work with Point shapefiles (Python GDAL Library)

texts:
def main1(shapefile):

    #making the shapefile as an object.
    input_shp = ogr.Open(shapefile)

    #getting layer information of shapefile.
    shp_layer = input_shp.GetLayer()

    #pixel_size determines the size of the new raster.
    #pixel_size is proportional to size of shapefile.
    pixel_size = 1

    #get extent values to set size of output raster.
    x_min, x_max, y_min, y_max = shp_layer.GetExtent()

    #calculate size/resolution of the raster.
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)

    #get GeoTiff driver by
    image_type = 'GTiff'
    driver = gdal.GetDriverByName(image_type)

    #passing the filename, x and y direction resolution, no. of bands, new raster.
    new_raster = driver.Create(output_raster, x_res, y_res, 1, gdal.GDT_Byte)

    #transforms between pixel raster space to projection coordinate space.
    new_raster.SetGeoTransform((x_min, pixel_size, 0, y_min, 0, pixel_size))

    #get required raster band.
    band = new_raster.GetRasterBand(1)

    #assign no data value to empty cells.
    no_data_value = -9999
    band.SetNoDataValue(no_data_value)
    band.FlushCache()

    #main conversion method
    gdal.RasterizeLayer(new_raster, [1], shp_layer, burn_values=[1])

    #adding a spatial reference
    new_rasterSRS = osr.SpatialReference()
    new_rasterSRS.ImportFromEPSG(2975)
    new_raster.SetProjection(new_rasterSRS.ExportToWkt())
    return gdal.Open(output_raster)

This is my code adapted from the Python Geospatial Analysis Cookbook.
I'm trying to convert a vector .shp to a raster .tif. The code runs and works for line and polygon vector files, but for point shapefiles, it creates an empty black tif file.

Here is the input:

This is the output

So what changes do I make in the code so that it works for point shapefiles too?


texts: ArcGIS compute nearest network neighbour for each point

texts:
I have a set of points and I would like to calculate a value for each point to indicate how 'clustered' to other points that point is.

Distance to nearest neighbour (through a network) would do fine.  Other clustering measures with significance would be good too: but Getis Ord GI and Local Moran's I rely on the points having attribute values, which mine don't.

How can you do this in ArcGIS, produce a measure for each point that tells you something about how clustered that point is, without attribute values?


texts: Using ArcGIS Javascript API and need help understanding asynchronous programming in accessing the fields property of a FeatureLayer Object

texts:
I'm creating a new FeatureLayer, and I want to populate the popup with a list of all the fields. Here is part of my code:

var featureLayer = new FeatureLayer(url, {
    mode: FeatureLayer.MODE_ONDEMAND,
    outFields:["*"],
    infoTemplate: popupTemplate
});
var fieldArray1 = featureLayer.fields;
console.log(fieldArray1); // output: undefined

// But, this gives me what I want:
setTimeout(function(){
    var fieldArray2 = featureLayer.fields;
    console.log(fieldArray2); //output: an array of Field Objects
},5000);

I believe one of my problems is that I don't fully how to handle the fact that javascript is asynchronous. I realize that var fieldArray1 = featureLayer.fields; is executing before the FeatureLayer Object finished being created.

What is the correct and best way to code what I want to do, which is - "after the new FeatureLayer Object is created, set fieldArray1 = featureLayer.fields?"


texts:
In the July 2017 edition of Swiss Review, a publication for Swiss living abroad, the cover feature was on Swiss Cartography.

With my interest in GIS this obviously sent me straight to the article.

See http://www.revue.ch/fileadmin/revue/Ausgaben/2017/04/SRV_1704_EN.pdf page 6.

One of the sections looks at how the usual Swiss, somewhat humourless approach to attention to detail and precision, has lead to some "skilfully crafted jokes" being introduced into the official published maps.

Here are those ones presented that escaped the "sharp eyes of proofreaders"

I am putting together a training presentation for QGIS and these sorts of jokes make for nice light-hearted, ice-breaker content that helps the flow of what can be a technical presentation.

Does anyone have examples of cartographic jokes that they have come across, or even better, introduced themselves into official, and not-so official published maps that can help liven up a QGIS presentation?


texts:
I am looking for a method to export data (Feature Layers) from my Company's ArcGIS Online account.  I see that there is a option to export data in various formats in the Gallery section, but we want to make this export as seamless as possible.  The objective is to input this data into Pointerra to view GIS data that is stored in ArcGIS Online.

How can I make my ArcGIS Online data available via WMS so a software like Pointerra can read my data?

Pointerra can also stream in TMS, WMTS, ESI REST Services (Not primary option).


texts: How to concatenate strings separated by dot in PyQGIS?

texts:
I would like to concatenate the string of field1 with field2 separated by dot and so on. However, when I put the "." the attribute values remain null. How to separate string by dot?

e = QgsExpression('concat((lpad("field1",2,0)),".",(lpad("field2",2,0)),".",(lpad("field3",3,0)),".",(lpad("field4",4,0)))')


texts:
I am using ArcGIS Desktop.

I have 33 DEM tiles that I am trying to mosaic together using this tool,

but I keep getting strange results.

All of the DEM pixel depths/sizes are 32 bit float, 1 band and are in NAD83, which I have been specifying in the tool. I'm saving it to a gdb and not specifying a cellsize.  I am also leaving the mosaic operator and colormap modes as the defaults.  This is the result I am getting:

Wondering if anyone knows why this is happening.

I tried using a mosaic dataset, but I am having issues with performing additional manipulation of it after I added all of the tiles.  I would prefer to use the mosaic to new raster as I know the additional analysis and manipulation I need to perform has worked in the past.


texts: Question on how to identify (incorrectly tagged) bicycle routes from OSM

texts:
My workmates and I are trying to use OSM data about bicycle routes and other touristic infrastructure for the assessment of cultural ecosystem services. First, we compared the OSM data to official data to get an idea of completeness and correctness of the data. However, the completeness of the data for our region near the Bavarian-Czech border is surprisingly low, many cycleways seem to be missing in OSM.

Interestingly, there are some roads and tracks in OSM, which might be used by cyclists but they don´t have any bicycle tags. Some are tagged as "hiking" for example. I already read the OSM guidance on how to tag new data but at least some of these routes are obviously not correctly tagged.

Can you please help me and explain to me how others extract cycle routes from OSM data without missing the incorrectly tagged parts? It must be possible because cycling apps, which are based on OSM, don´t have as many gaps. Are they corrected manually?

I hope, you can help me.


texts:
For a project I want to create a little program that does simple things like export a map from the python command line, but also users are able to edit said map (without having to open arcmap) e.g change font text, replace a street name with another street name, change map layout. I know you can do something similiar in arcpy.mapping from the command line - but what if I use gui tools to create a small program that allows users to do so by pressing buttons (e.g one button exports map, I could then upload a few example maps that can be worked on) - I'm looking at using tkinter.

Does this sound overly complex?

I can remove the gui tools part if its too hard - is it possible to test out python programs somewhere?

like can I upload it somewhere and test it out and do I need html knowledge?
on a second note - does anyone know a way to replace a graphic in all maps - for example, if I want to replace the company logo of all my maps with another logo I've made - I've googled and can't seem to find an answer to it.


texts:
I'm following this tutorial. I installed the “GDAL Complete” Framework from kyngchaos.com and I tested the installation on Mac. I have GDAL 1.11.2, released 2015/02/10.

When I am trying to run this command:

ogr2ogr -f MySQL MySQL:dbname,user=root,password=password \
      -nln world_borders -nlt MULTIPOLYGON \
      -update -overwrite \
      -lco ENGINE=MyISAM \
      -lco MYSQL_FID=ogr_fid -lco GEOMETRY_NAME=geometry \
      /path/to/TM_WORLD_BORDERS-0.3.shp

With my settings it would be:

ogr2ogr -f MySQL MySQL:mapas,user=xxxx,password=xxxx \
      -nln test -nlt MULTIPOLYGON \
      -update -overwrite \
      -lco ENGINE=MyISAM \
      -lco MYSQL_FID=ogr_fid -lco GEOMETRY_NAME=geometry \
      /Users/Seph/Downloads/TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp

I got this error:

Unable to find driver `MySQL'.
The following drivers are available:
-> `ESRI Shapefile'
-> `MapInfo File'
-> `UK .NTF'
-> `SDTS'
-> `TIGER'
-> `S57'
-> `DGN'
-> `VRT'
-> `REC'

I have created the table on my database before running the command.

UPDATE
I installed Homebrew I and ran this command:

brew install gdal --with-mysql

Now I have the supported format:  -> "MySQL" (read/write)

But now I have this error:

Macbook-pro:~ Seph$ ogr2ogr -f MySQL MySQL:mapas,user=xxxx,password=xxxx\
      -nln test -nlt MULTIPOLYGON \
      -update -overwrite \
      -lco ENGINE=MyISAM \
      -lco MYSQL_FID=ogr_fid -lco GEOMETRY_NAME=geometry \
       /Users/Seph/Downloads/TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp
      ERROR 1: MySQL connect failed for: mapas,user=xxxx,password=xxxx
      Can't connect to local MySQL server through socket '/tmp/mysql.sock' (2)
      ERROR 1: MySQL connect failed for: mapas,user=xxxx,password=xxxx
      Can't connect to local MySQL server through socket '/tmp/mysql.sock' (2)
      ERROR 1: MySQL connect failed for: mapas,user=xxxx,password=xxxx
      Can't connect to local MySQL server through socket '/tmp/mysql.sock' (2)
      ERROR 1: MySQL driver doesn't currently support database creation.
      Please create database before using.
      MySQL driver failed to create MySQL:mapas,user=xxxx,password=xxxx

I am using phpMyAdmin from Xammp


texts: Using NDVI values as covariate in survival analyses

texts:
I am trying to become more familiar with processing NDVI values in GIS. However, despite its usefulness, it seems to have some limitations. Namely, creating NDVI values across a large scale can be problematic if shading (i.e. cloud cover) and other lighting conditions are not constant especially when comparing NDVI values of the same area over time (if I am understanding this correctly). I have local digital images of 1 meter square plots taken from a high resolution hand-held camera (and with a hand-held NDVI camera) of vegetation plots at a height of ~5 feet. I would like to generate an average NDVI value for these images and use them as covariates in a survivial model of my target species. Ideally, I would upload these photos as TIFFs into ArcMap for processing. I have seen people use the red and infra-red bands from Landsat images to generate an NDVI value for the raster layer. Could I apply this to my photos?

Also, is it acceptable to use the NDVI value at each location as an independent predictor variable? Or, will the local NDVI values be confounded by variation in lighting conditions even if I'm not comparing the values from each image directly? I standardized the area photographed, time of day (10am-2pm) photos were taken, and camera settings.


texts: Extracting Road Terrain Properties using DEM / LIDAR

texts:
I want to extract following road terrain data from DEM but I'm not sure if its possible. Can i extract these attributes using LIDAR if not DEM ?

These are the data fields I need:

Horizontal curve location, Horizontal curve length,    Grade in direction of travel,   Roadway elevation,  Roadway elevation,  Vertical curve location,    Vertical curve length,  Available sight distance,   Road cross slope, Shoulder cross slope

Any information would be helpful.


texts: How do I find the intersection of two polylines in PostGIS?

texts:
I have a set of road conditions reported for small segments of road, and I am hoping to determine the specific subset of each segment which intersect with a route along the road so that I can draw the intersection of the route and road conditions on a map.

Is there a way to determine the specific points of intersection of two polylines in PostGIS in a way suitable for plotting on a Google Map?

I am new to PostGIS so please accept my apologize if I am getting any terminology wrong.


texts: Want a transformation box after a user draws a polygon

texts:
I have built an application which allows users to draw a polygon inside another polygon. The two polygons are in different layers for some application specific reasons. After the user draws a polygon, layer.events.register function is invoked to check if the polygon lies entirely inside the larger polygon. Now I want to modify the application and I have the following questions:

Q1) How can I add a transformation box to the polygon that is drawn by the user? I want the user to be able to rotate or scale the polygon after it is drawn.

Q2) How can I check that after the transformation(rotation or scaling) the polygon will still be inside the larger polygon?

** The polygon drawn by the user is in fact a rectangle inside another rectangle.

I have pasted some code below. Note that the layer named "BaseLayer" is where the user generated rectangles are added.

var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
var baseLayer = new OpenLayers.Layer.Vector("BaseLayer", {
            styleMap: new OpenLayers.StyleMap({
                "transform": new OpenLayers.Style({
                    display:"${getDisplay}",
                    cursor: "${role}",
                    pointRadius: 3,
                    fillColor: "white",
                    fillOpacity: 1,
                    strokeColor: "black"
                    },{
                 context:{
                    getDisplay: function(feature){
                        return feature.attributes.role === "se-resive" ? "none": "";
                        }
                    }
                })
            }),
            renderers: renderer
});

var transformControl = new OpenLayers.Control.TransformFeature(baseLayer, {
                renderIntent: "transform",
                rotationHandleSymbolizer: "rotate"
});
var drawPolygonFeature = new OpenLayers.Control.DrawFeature(baseLayer, OpenLayers.Handler.RegularPolygon,
            {'displayClass':'olControlDrawFeaturePolygon',
                 handlerOptions:{sides: 4,
                 irregular:true,
                 style: polygonStyle,
                 id: "processingPolygon"
                 }
            });
baseLayer.events.register("featureadded", feature, function(evt){
    // DOES SOME ERROR CHECKING HERE.
   //  transformControl.setFeature(evt.feature);
});

What I've tried so far:
I tried adding the transformation box to the feature inside the events.register() function (second last line in the code above) using the control.setFeature function, however this did not work. Adding a transformation box also invokes the events.register function which adds another transformation box and as a result invoking the events.register() again. This results in an infinite loop!

Edit:
Here is an example of the type of transformation box I want: http://openlayers.org/dev/examples/transform-feature.html


texts: I'm looking for an alternative to 'snakes method' for smoothing very big polygons, which I created from SRTM elevation datasets. The problem is that the smoothing is getting stuck after a while with default parameters. Can I change the parameters for better performance or should I choose another method for this?


texts: How can I use a tif and a SID to make a basemap in QGIS?

texts: I am trying to add a basemap in QGIS, projected correctly for PA.  I have been directed to PASDA.  When I click on the FTP download link, it gives me the option to download TIF data or SID data.  What do I want, and how will I add it to QGIS as a basemap?


texts: Investigating PostGIS functions to see what they are doing?

texts: Let's say you have a function like ST_Distance and I want to know how it works. The docs don't say much at all, how do I go about debugging this or figuring out how it ticks.


texts: How to manually add texts on the map in QGIS?

texts: Is there another way to put river, towns etc names beside labeling in QGIS? Can you manually type text determine size and place it where and as you wish on the map?


texts: Problem with CRS and proj.db using Reproject from Rasterio

texts:
I saw similar problems but I still can't manage with mine.

Some technical information:

I've installed rasterio version via anaconda (conda install rasterio): rasterio-1.2.1-cp37-cp37m-win_amd64.whl
and GDAL version: GDAL-3.2.2-cp37-cp37m-win_amd64.whl
I set environmental variable PROJ_LIB to C:\Users\user\anaconda3\Library\share\proj

Code

I try to use reproject function from rasterio.
with rasterio.open("C:\\Users\\user\\Desktop\\Sentinel_2\\S2A_OPER_MSI_L1C_TL_EPA__20160815T031336_A000462_T34UEV_B02.jp2") as scl:
B02=scl.read()
reproject(
    B02, tmparr,
    src_transform = scl.transform,
    dst_transform = aff,
    src_crs = scl.crs,
    dst_crs =scl.crs,
    resampling = Resampling.bilinear)
B02 = tmparr
print(B02.shape)

And I received error:
Missing src_crs
Then I checked my image EPSG code and paste it to code:
with rasterio.open("C:\\Users\\user\\Desktop\\Sentinel_2\\S2A_OPER_MSI_L1C_TL_EPA__20160815T031336_A000462_T34UEV_B02.jp2") as scl:
B02=scl.read()

src_crs = {'init': 'EPSG:32634'}

reproject(
    B02, tmparr,
    src_transform = scl.transform,
    dst_transform = aff,
    src_crs = src_crs,
    dst_crs = src_crs,
    resampling = Resampling.bilinear)
B02 = tmparr
print(B02.shape)

And I get error:  The EPSG code is unknown. PROJ: proj_create_from_database: Cannot find proj.db
I would be grateful for all the tips


texts:
I have a CSV file. I would like have this file inside GeoServer. I used "add postgis layers" in QGIS and I added the connection that I have in  PgAdmin III. I had the message that the connection is succesful.

I added the CSV file using "DBManager" in QGIS. The problem is that the section "public" inside "DBManager" is red (means error) and I don't have the CSV in PgAdmin III.

This is the error that I had


texts: Display raster with earthpy shows one color raster without values

texts:
I'm new to Earthpy, I am trying to display raster and on top of it I want to display polygon. The problem is that when I display the image it looks like blue rectangle and doesn't have different colors according to the NDVI values that the image has.
This is the original image:

and this is how it looks when I try to display both layers together:

This is my script:
#load shapefile
polygons=geopandas.read_file(r'very\nice\path\polygons.shp')

#change the projection of the shape to be the same as the raster one (I know the projection)
polygons.crs = {'init' :'epsg:4326'}
polygons=polygons.to_crs({'init': 'epsg:32721'})

#Load image
img=rasterio.open(r'My\nice\path\Raster.tif')
extent = rio.plot.plotting_extent(img)

#Have only one band to display
imageband=img.read(1)

#display raster and on top of it the polygons:

fig, ax = plt.subplots(figsize=(10, 8))

ep.plot_bands(imageband, cmap='terrain',
              extent=extent,
              ax=ax,
              title="Raster Layer with Shapefile Overlayed",
              cbar=False)

polygons.plot(ax=ax, alpha=.8)

ax.set_axis_off()

My end goal is that when I plot the polygons, the background image (e.g raster) will look the same as it looks when I upload it at the beginning
I'm also open to other solutions using other libraries


texts: ArcCatalog and python crash when running ConvertTimeField_management

texts:
I have a script which reclassifies a tif file, adds a new field to the tiff, and is supposed to parse the date from the filename and add it to the new date field.  Everything works as it should, except parsing the date.  When the script runs this line, python crashes, either runnig the script through ArcCatalog or python command line:

arcpy.ConvertTimeField_management("C:\Users\OuelletteMS\Desktop\Ice_Data_from_JohnIacozza\1978\wa_19780629_reclass.tif", 19780629, "yyyyMMdd", "Date")

This example is using variable values as they are assigned by the script (print the variables, copy/paste and add quotes).  The full script is below.

I get no error messages, and commenting out this one line shows that everything else works.  I got this line as a snippet after a successful run of the tool in ArcMap.
What am I missing?

I'm using ArcGIS 10.1 SP1 and Python 2.7

import arcpy, os, sys
from arcpy import env
arcpy.CheckOutExtension("Spatial")

# Allow for file overwrite
arcpy.env.overwriteOutput = True

# Set the workspace directory
path = "C:\Users\OuelletteMS\Desktop\Ice_Data_from_JohnIacozza"
env.workspace = path

# Get the list of the featureclasses in the workspace directory
listFCs = arcpy.ListFiles("wa_197806*.tif")

# Loop through each file and perform the processing

for fc in listFCs:
    print str("processing " + fc)       # displays the file that is currently being handled

    shpname = str(fc)                   # converts filename to a string
    yearVal = shpname[3:7]              # year processed
    newName = shpname[0:11] + "_reclass.tif"    # name for reclassified raster
    newPath = path + "\\" + yearVal     # sets the path for the reclassified file
    outfile = newPath + "\\" + newName  # specifies reclassified raster path
    vf = "Value"                        # sets value field to be used
    strDate = shpname[3:11]             # extracts characters 3 through 11 from the file name

# Check to see if a directory for that year exists.  If not, make one
    try:
        os.makedirs(newPath)
    except OSError:
        if not os.path.isdir(newPath):
            raise

# Reclassify raster
    arcpy.gp.Reclassify_sa(fc,vf,"1 5 1;6 99 0;NODATA 0",outfile,"NODATA")

# Adds and populates date field
    field2 = "Date"                     # Field for storing the date as a date
    arcpy.AddField_management(outfile,field2,"DATE")
#    arcpy.ConvertTimeField_management(outfile, strDate, "yyyyMMdd", field2)

# Cleanup
del listFCs, fc, shpname, yearVal, newName, newPath, outfile, vf, strDate, field2

#
arcpy.CheckInExtension("Spatial")
print "done"


texts: Using legend with dijit/form/DropDownButton in ArcGIS API for JavaScript?

texts:
I am very new to the ArcGIS API for Javascript and am by no means a programmer. Nonetheless, I have been fairly successful at finding existing code and modifying it to my needs. That is, until now.

I want something that I think should be simple and straightforward. However, despite looking high and low, I cannot find a good code sample. I would like a button that, when clicked, shows or hides a legend (or TOC). Indeed, I have seen multiple examples of this in WebAppBuilder applications and elsewhere, but it has not been easy to tease through the code. Surely, some others have done this.

It seems to me that my answer should be in "dijit/form/DropDownButton." However, this newbie cannot seem to get it to work. My guess is I am overlooking the obvious (or, perhaps, I am approaching it all wrong).

I have attached a snippet of the code sample below.

 require([
    "esri/map",
    "esri/layers/FeatureLayer",
    "esri/InfoTemplate",
    "esri/renderers/UniqueValueRenderer",
    "esri/symbols/SimpleMarkerSymbol",
    "esri/Color",
    "esri/dijit/Legend",
    "dijit/form/DropDownButton",
    "dojo/_base/array",
    "dojo/domReady!"
  ], function(Map, FeatureLayer, InfoTemplate, UniqueValueRenderer, SimpleMarkerSymbol, Color, Legend, DropDownButton, arrayUtils) {
    map = new Map("map", {
    basemap:"gray",
    center: [-98.435731, 35.222876],
    zoom: 7,
    logo: false,
    slider: false
  });

    var layer = new FeatureLayer("//services.arcgis.com/V6ZHFr6zdgNZuVG0/ArcGIS/rest/services/NationalParkStats2013/FeatureServer/0");

    var renderer = new UniqueValueRenderer(null, "Type");

    var symbol1 = new SimpleMarkerSymbol();
    symbol1.setColor(new Color("#ed5151"));
    renderer.addValue("National Park", symbol1);

    var symbol2 = new SimpleMarkerSymbol();
    symbol2.setColor(new Color("#149ece"));
    renderer.addValue("National Monument", symbol2);
    layer.setRenderer(renderer);

    map.addLayers([layer]);

        //add the legend
      map.on("layers-add-result", function (evt) {
        var layerInfo = arrayUtils.map(evt.layers, function (layer, index) {
          return {layer:layer.layer, title:layer.layer.name};
        });
        if (layerInfo.length > 0) {
          var legend = new Legend({
            map: map,
            layerInfos: layerInfo
          }, "legend");
          legend.startup();
        }

      });

  });

My HTML body then has these divs:

It certainly would be nice if this worked!

No really

        Legend


texts: Select Layer By Attribute Python Script Tool?

texts:
I'm trying to script another Select by Attributes tool. I want to be able to find a specific building using a field and zoom in to the location. However, my tool needs to allow a user to input any field or any value. So far this is what I got.

#Import Modules
import arcpy

arcpy.env.workspace = "C:\Users\pierrej\Desktop\GIS Data"
arcpy.env.overwriteOutput = True

#Set to current mxd and dataframe
mxd = arcpy.mapping.MapDocument ('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]

try:
    # Set the tool parameters
    InputFeatureClass = arcpy.GetParameterAsText(0)
    InputField = arcpy.GetParameterAsText(1)
    InputValue = arcpy.GetParameterAsText(2)

    # SQL expression used in selecting a feature
    #Create the query
    qry = "\"InputField\" = InputValue"

    arcpy.MakeFeatureLayer_management(InputFeatureClass, "Layer")

    #Call the Select Layer by Attribute tool
    arcpy.SelectLayerByAttribute_management("Layer", "NEW_SELECTION", qry)

    #Zooming to a selection set for the specific layer
    df.zoomToSelectedFeatures()
    df.scale = 2500000
    arcpy.RefreshActiveView()

    # Report a success message
    arcpy.AddMessage("All done!")

#Print an error message in the event of a problem:
    except:
    print "An error occurred during selection"

My parameters are the following:

InputFeatureClass = Data Type: Feature Layer, Required, Input, No Multivalue

InputField = Data Type: Field, Required, Input, No Multivalue, Obtained from Input Feature Class

InputValue = DataType: String, Required, No Multivalue

The script runs and I don't have any error message but the tool doesn't select any thing. How can I fix this?


texts: Replicating Conditional popup of ArcGIS Online in QGIS2WEB?

texts:
I stumbled across this ESRI post about Conditional Field display with Arcade in Pop Ups. I am not well versed in HTML, but I could follow along well enough to understand the post.

However, how would I implement this in QGIS, specifically through the QGIS2WEB plugin?

I am examining the html output of the QGIS2WEB plugin (for both Open Layers and Leaflet), and an obvious intervention point is not jumping out at me.

Popups can sometimes have empty fields when there is no data, and I would like to hide these fields when there is no data!


texts: Shapes of this entity type not allowed in this layer when CreateFeature() is called

texts:
I have the following code that calls on a remote geodatabase server. Everything is connected and working, but when I try to add a point and use featureclass.CreateFeature();, it shuts down.

At first I realized it was esriMultiPoint, but I fixed this, I remade the command lines, and created the shp2sde again. And this time I can confirm the server says it is esriGeometryPoint.

IWorkspaceEdit wse = ws as IWorkspaceEdit;
wse.StartEditing(true);
wse.StartEditOperation();
Log.Write("Feature class=" + fClass.FeatureType.ToString() + " Fieldtype=" +fClass.ShapeType.ToString());
//Get a reference to the new point object
ESRI.ArcGIS.Geometry.IPoint com_point = (ESRI.ArcGIS.Geometry.IPoint)(serverContext.CreateObject("esriGeometry.Point"));
Log.Write("Create esriGeometryPoint");
((ESRI.ArcGIS.Geometry.IGeometry)com_point).SpatialReference = ((IGeoDataset)fClass).SpatialReference;
IGeometry geometry_point = com_point as IGeometry;
Log.Write("Check geometry point type=" + geometry_point.GeometryType);
//input coords
object tempObj;
paraMap.TryGetValue("Longi", out tempObj);
com_point.X = double.Parse(tempObj.ToString());
paraMap.TryGetValue("Lati", out tempObj);
com_point.Y = double.Parse(tempObj.ToString());
Log.Write("about to create feature");
IFeature f = fClass.CreateFeature();
Log.Write("Create feature"); //CODE NEVER GETS HERE
f.Shape = com_point;

Here's my commandline to create the table in the database... (I put dots for excluding some info).

shp2sde -o create -l mynewdata3,shape -f C:\..bla.shp -a all -i sde:sqlserver:..... -C FID,SDE,1 -s ...... -D ...... -e p -u ..... -p ......

EDIT, should have used:

shp2sde -o create -l mynewdata3,shape -f C:..bla.shp -a all -i sde:sqlserver:..... -C FID,SDE,1 -s ...... -D ...... -e pn -u ..... -p ......

The log output and ERROR:

Feature class=esriFTSimple Fieldtype=esriGeometryPoint
Create esriGeometryPoint
Check geometry point type=esriGeometryPoint
about to create feature
at ESRI.ArcGIS.Geodatabase.IFeatureClass.CreateFeature()\r\n Shapes of this entity type not allowed in this layer

EDIT: I checked with logs and saw that the feature class being used, is correct feature class too.

EDIT 2: Solved

Details about the error:
http://support.esri.com/en/knowledgebase/techarticles/detail/16135

Solution:
http://forums.esri.com/Thread.asp?c=159&f=1707&t=232881


texts: Retrieving MXD Display view page units with ArcPy?

texts:
I've read some articles about translating map units into page units for getting distances in an MXD.  I have not seen a post about how to get the page units programmatically from ArcPy.  Is it possible to in any version of ArcMap 10.x to get the page units programmatically?  DisplayUnits and mapUnits from dataframe do not suffice.

Edit:
 This post   only answers how to translate units between map and page units.  It does not answer how to retrieve with arcpy page units (cm, in, etc.).
I am not looking for the answer to how to translate units from map distance to page distance and back.  I am simply asking about how to retrieve page units.


texts: `change:center` event not fired on drag animation

texts:
I'm quite new in OpenLayers, and I want to act when the map is moved.

So simply listening to the change:center event of the map's view:

map.getView().on('change:center', function(evt){
   console.log(map.getView().getCenter());
}

But I've noticed the event is not fired after we finish the map panning with the mouse - and the map is still moving and slowing down until it stops completely.

Is there another way to know when the map extent is changed?


texts: Generate an sql WHERE clause from selected features


texts:
I'm trying my first GIS project with cartopy and matplotlib and I'm having some trouble displaying some shapefiles.

The shapefiles were downloaded from here: http://data.vancouver.ca/datacatalogue/localAreaBoundary.htm

The code I'm trying to use has worked for other shapefiles but not for this set.

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

fname = "local_area_boundary.shp"
adm1_shapes = list(shpreader.Reader(fname).geometries())

ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_geometries(adm1_shapes, ccrs.PlateCarree(),
                  edgecolor='black', facecolor='gray', alpha=0.5)
east = -123.0
west = -123.3
north = 49.1
south = 49.4
ax.set_extent([west, east, south, north], ccrs.PlateCarree())
plt.show()

Again, I'm 100% new to GIS and cartopy (but fairly experienced in Python) and I'm sure I'm missing something easy. Any help debugging this would be great.

Edit: I've modified the code so the shapefile is displayed, but it now has the wrong lat/lon coordinates

fname = 'local_area_boundary.shp'
adm1_shapes = list(shpreader.Reader(fname).geometries())

f, ax = plt.subplots()
ax = plt.axes(projection=ccrs.Mercator())

ax.add_geometries(adm1_shapes, ccrs.Mercator(),
                  edgecolor='black', facecolor='white', alpha=0.5)
a = 483644.82757398
b = 5449579.62327763
c = 498313.01571369
d = 5460349.22564139
ax.set_extent([a,c,b,d ], ccrs.Mercator())

I would expect the lat/lon to be centred around 49.2N/-123.1E, but instead it's 44.1N/4.3E. Can you tell if this is just bad data, or if I'm still doing the transformation wrong?


texts: Is there an easier way to edit polygons in ArcMap?

texts: I am updating woodlot polygons in ArcMap 10. Instead of editing individual verticies, is there a way to edit that is more like drawing part of a new shape (say if the woodlot has expanded only on one side). I'm looking for something similar to the trace tool but I'm working off a tiff image so trace doesn't work. I just want to trace the woodlot and have the existing polygon update to that new shape. Anything to make this task faster would be much appreciated!


texts:
I´m making a GIS-Portal with GeoNode 2.6 (running in ubuntu 16.04 server) and I want to change the look of the home page following this instructions: http://docs.geonode.org/en/master/tutorials/admin/customize_lookfeel/customize/theme_admin.html (I´m a beginner at this). I have copied the "site_index.html" and "site_base.html" from
/usr/local/lib/python2.7/dist-packages/geonode/contrib/geosites/site_template/templates/
to
/etc/geonode/templates/
I edited both html files.
Then I run the collectstatic command and restart apache2, but i see no changes in the homepage.

in local_settings.py is set:

TEMPLATE_DIRS = (
    '/etc/geonode/templates',
    os.path.join(GEONODE_ROOT, 'templates')

Any ideas of what i´m missing to make it work?

I also tried editing the index.html file without any luck. It only changes when editing /usr/local/lib/python2.7/dist-packages/geonode/templates/index.html


texts: Looking For Free Global Bedrock Geology Data

texts: I've searched around for a freely available global bedrock geology data set, but I have yet to find one.  Anyone come across such dataset?


texts: ArcGIS Python Tool - Importing Custom Script Into ToolValidator Class

texts:
I had posted a question last week about customizing a ToolValidator class and got some very good answers.  In working with the proposed solutions, I have created a custom module that performs queries on a db, and will be called by both the ToolValidator class (to provide values for the drop-down lists) and later in the geoprocessing script (to get other parameters based on items selected in drop-down lists).  However, I can't seem to actually call the custom module in the ToolValidator class.  I have been trying to append to path with no luck.  When I try to apply these changes to the script, I get a runtime error:  [Errno 9] Bad file descriptor.  If I comment out import line, no errors.

sys.path.append('my_custom_module_directory')
import my_custom_module

Many of you may be asking why don't I just implement a custom tool with ArcObjects.  The reason is that my end users don't have the privleges necessary to register ANY dlls on their computer.

UPDATE:  This was happening to me in ArcGIS 10.  Interestingly enough, I was initially appending to the path inside of the initialiazeParameters function of the ToolValidator class.  If I do the append outside (i.e., on top of) of the ToolValidator class, everything works as expected.

sys.path.append('C:/Working/SomeFolder')
import somescript -------->THIS WORKS

class ToolValidator:
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    import arcpy
    sys.path.append('C:/Working/SomeFolder')
    import somescript -------> THIS DOESNT WORK
    self.params = arcpy.GetParameterInfo()

UPDATE 2:  I think I found the true cause of my problem.  In the code snippets in this post, I have been appending what appear to be real paths (ex C:/Working/SomeFolder) to the sys.path.  In my actual ToolValidator class, I was constructing a relative path using os.path.dirname(__file__) +"\my_special_folder...".  I was anticipating that os.path.dirname(__file__) would return the path of the toolbox, since that contains the ToolValidator class.  I have come to find this is not the case.  As far as I can tell, the ToolValidator class is never actually written to a .py file, and I speculate the the this code is passed to the python interpreter in memory, so __file__ is useless, or some temp script is persisted and then execfile(path_to_script) is called, again rendering __file__ useless.  I'm sure there are other reasons I am missing as well.

Long story short, if I use a hardcoded path, the sys.append works anywhere, relative paths don't work so well in the ToolValidator class.


texts: How to get current application culture using ArcObjects?

texts:
Without looking at the registry, is there a way to get the current application culture in C#?

I have a Add-in and want to get the language from the ArcMap hook using ArcObjects 10.

(I'm using C# to develop the addin.)


texts: Copying features from if field value is using ModelBuilder

texts:
I am currently trying to understand ModelBuilder and have decided to break apart a big project into smaller ones to understand how everything works. Currently I am stuck on trying to copy features from an if field value is logical tool.

The problem I am having is that I cant seem to copy the features from my true or false values. I've tried using copy features and but it wants an input and does not take the true or false is one. If I just in the point features I am using it will just spit them back out since its just copying those fields.

Does anyone have any ideas on what other tools to use or an input to create?


texts: Changing geographic coordinate system of feature class in ArcMap?

texts:
How do I make these two feature classes compatible in ArcMap 10.4?

(1)  -->
Data Type:  Shapefile Feature Class
Shapefile:  C:\Users\craigugoretz\Desktop\CA_counties\CA_counties.shp
Geometry Type:  Polygon
Coordinates have Z values:  No
Coordinates have measures:  No

Geographic Coordinate System:   GCS_WGS_1984
Datum:  D_WGS_1984
Prime Meridian:     Greenwich
Angular Unit:   Degree

(2) -->
Data Type:  Shapefile Feature Class
Shapefile:  C:\Users\craigugoretz\Desktop\NHD\WBDHU12.shp
Geometry Type:  Polygon
Coordinates have Z values:  No
Coordinates have measures:  No

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree


texts: Automatically fixing multiple errors using Fix topology error tool?

texts:
I am an ArcGIS 10 user.

I want to solve multiple errors using the fix topology error. I select all my features, I click right but I don't have  the correcting options (Snap, trim, extent), I only have Delete.

Does anyone know why?


texts: Is Tiger/Line 2010 Data Inaccurate? What's more reliable?

texts:
I followed the instructions at:
http://www.letseehere.com/postgis-geocoder-using-tiger-2010-data

To get postgres working with postgis, but some of the addresses are off 500 ft to as much as 1 mile in my random sampling. Do you know of any more reliable data sets?


texts:
I've got some old maps in dxf format. QGIS imports them correctly but with a nonsense scale. It can export the data it gets to shp of course. Does anyone know how to rescale vector data so it can be related to other maps? It doesn't have to be done in QGIS as long as the method is open source and it works on a format QGIS can handle (SHP is best, of course).

Note that there seem to be two scales in SHP: the data I hand-imported using WKT shows the correct scale on the mouse, but the map scale indicated is nonsense, and print composer is affected the same way. I'd like to correct both types of problem.


texts: How to split a GPX track file into several files of N trackpoints each?

texts:
I have some GPX track files, each containing tracks with thousands of trackpoints.

I need to split those files into multiple ones, each with 500 (or more, accordingly) points each.

Unfortunately I cannot simplify those tracks, I really need to split them.

As long as I could find, gpsbabel isn't able to do so. I am open to use any tool either on Linux (preferably) or Windows as well.


texts: Finding Mean Distance to 3 Locations from POI using ArcGIS Network Analyst

texts:
I am wondering how to complete the following task using Network Analyst in ArcMap 10.7.1.
I have a point file of fire station locations and also a point file of office locations. I would like to find the mean distance to 3 fire stations from all of the office buildings. I was thinking of doing Closest Facility but I think I would need to do that for all the fire stations which I can't do because there's too many.
Any suggestions?


texts: Selecting text in parentheses in field using ArcMap

texts:
I have imported CAD annotation and changed it to point feature, and I have the "text string" field in table attribute. Now I want to filter any records that contain the text in parentheses.
For example: the text is "100PVC (2005)" => I want to get "2005"
But the field is more complicated than this example, because some records just like "(05-07)".
I have tried the index function with field calculator but it did not work. I also could not select records which contain the parentheses by using the Select by Attributes. My expression was:

[TextString] LIKE '%(____)%'

When I verify that it was successfully but no records were returned.

How can I solve this problem?


texts: Making points count automatically when zooming out in QGIS?

texts:
Can I make the point count itself when there are too much stacked up points in one area when I zoom out my map?

I want my output e.g. Area of ganli points to be like this when i zoom out

Water point functionality image is not mine i just post it for sample.


texts: Multilinestring to Linestring PostGIS 2.3

texts:
I am trying, as the title suggests, to transform a geometry that was a MultiLinestring into a Linestring. I first create a copy of the table I want to change, as a temporary table from which I can UPDATE the geom of the real one after altering the geometry. I then use this following query:

update cablage_pays_gex.cable set geom = NULL;
alter table cablage_pays_gex.cable alter column geom type geometry(Linestring, 2154);
update cablage_pays_gex.cable
set geom = st_linemerge(cable_temp.geom)
from cable_temp
where cable_temp.id_cable = cable.id_cable;

But I get this error:

Geometry type (GeometryCollection) does not match column type (LineString)

What does it mean? What do I do wrong?


texts: Merging multiple satellite bands for science data set rasters?

texts: I am using MODIS satellite imagery data for one specific derived band. The problem I am running into is that it is a MCD19 SDS layer, which contains different bands for each orbit overpass at a specific location (up to 15 times). Because of this, I was wondering if it was possible to merge multiple bands together into one band to display in stretch symbology and average any overlapping imagery. I am currently using ArcGIS Pro but not able to find any geoprocessing or raster tools that fit this criteria for the TIFF images I am working with.


texts: Loop to sample random points within several polygons in R

texts:
I have a shapefile with 6 polygons, each one with a specific id.

I want to create random points within these polygons, but with different sample size (each id has a specific value).

The file with the polygons has the following structure:

> HR95_all
Object of class "SpatialPolygonsDataFrame" (package sp):

Number of SpatialPolygons:  6

Variables measured:
   id      area
0  40 257879080
01 41  53136845
02 42  64198336
03 43 186475935
04 44  27405855
05 45  53485801

And the data.frame with the number of points per id is a simple data.frame:

> AV_locs
   id  num
40 40 1010
41 41 1222
42 42  493
43 43 1111
44 44  369
45 45  598

I would like to create a loop that could do this automatically, and I've tried this:

for (i in 1:length(unique(HR95_all$id))) {
             id <- unique(HR95_all$id)[i]
             Rnd <- spsample(HR95_all[[id]], AV_locs[which(AV_locs$id==id), "num"] , "random")
             writeOGR(Rnd,dsn=getwd(), layer=paste0("Rnd_", id), driver = "ESRI Shapefile",overwrite_layer=TRUE)
     }

But I obtain the following error:

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘spsample’ for signature ‘"factor"’

I've tried other options, like including n = AV_locs[which(AV_locs$id==id), "num"] (following the help page), but the error persists.

Any idea on how to overcome this?

----- UPDATE -----

I've tried a different loop, by subsetting the HR95_all before applying the spsample():

for (i in (unique(HR95_all$id))) {
        HR <- HR95_all[HR95_all$id==i,]
        Rnd_locs <- spsample(HR, n = AV_locs[which(AV_locs$id==i), "num"], "random")
        writeOGR(Rnd_locs,dsn=getwd(), layer=paste0("Random_", i), driver = "ESRI Shapefile",overwrite_layer=TRUE)
}

The previous error is gone but a new one showed up:

Error in writeOGR(Rnd_locs, dsn = getwd(), layer = paste0("Random_", i),  :
  obj must be a SpatialPointsDataFrame, SpatialLinesDataFrame or
    SpatialPolygonsDataFrame

I've tried to add a as(Rnd_locs, "SpatialPointsDataFrame") before, but it still doesn't work; I recieve this error:

Error in writeOGR(Rnd_locs, dsn = getwd(), layer = paste0("Random_", i),  :
  NA NULL unknown data type

Any thoughts?


texts: ArcGIS - Getting Python toolbox to function properly

texts:
as the header says, I do have a problem in getting my toolbox to function the way I want it to. In the following my parameter definition is shown, where from I would expect to get only feature classes of type point to be shown.

    def getParameterInfo(self):
        #Define parameter definitions

        # Input Features parameter
        in_workspace = arcpy.Parameter(
            displayName="Input Workspace",
            name="in_workspace",
            datatype="DEWorkspace",
            parameterType="Required",
            direction="Input")

        # FeatureClasses parameter
        in_fc = arcpy.Parameter(
            displayName="Input FC",
            name="classes",
            datatype="DEFeatureClass",
            parameterType="Required",
            direction="Input")
        in_fc.parameterDependencies = [in_workspace.name]
        in_fc.filter.list = ["POINT"]

        parameters = [in_workspace, in_fc]

        return parameters

However, the resulting dropdown list shows all kind of content, e.g. point, line, and polygon feature classes. Even feature datasets are listed...

I got the suggestion to use GPFeatureLayer instead of DEFeatureClass, but as far as I understand Esri's resource center it looks to me as if GPFeatureLayer is in my case actually a detour to the datasets requested. When I tried GPFeatureLayer then there was no difference at all to the use of DEFeatureClass.

By my understanding, GPFeatureClass would be  the way to go, FeatureDatasets should not appear in the list, and the filter set for the list should only populate the list of FeatureClasses with datasets of type point.

Does anybody of you can help me with this issue?


texts: Reading .sde connection file without Esri softwares?

texts: I would like to be able to read connection properties from an sde file without any Esri dependencies such as Arcpy, ArcObjects etc. Is this possible?


texts: How to run MapServer on an existing XAMPP installation?

texts: How can I use MapServer on existing Apache installation? I don't want to install Apache server bundled with MapServer.


texts: Failed executing "Table to Excel" in ArcGIS

texts:
Can't export polygon shapefile into Excel sheet.
When executing Table to Excel in ArcGIS shows this error:

  ERROR 000576: Script associated with this tool does not exist.

  Failed to execute (TableToExcel).

What might be the cause of this?


texts: How to update my map with the restored features from JSON format?

texts:
I'm trying to implement a CRUD that have a map OpenLayers.

First, I drawed and saved the features with writeFeatures().

Second, I'm trying to show the polygon that I drawed, and here is my problem:

After I set the map.setLayer my map doesn't show the features.
I would like that when I select the product in the list, my map show the features that I read from my database, in format JSON. I saved using WriteFeatures();

Could I need a loop to iterate each item of array?

What am I doing wrong?
My configuration:

Laravel v6.18.22 with bootstrap;
OpenLayers v6.3.1

In my HTML has:
<div id="map" class="map" </div>

  <script>
    const getFeatur = @json(route('homeController.getLavouras'));
  </script>

 <script type="application/javascript" >
  my_map.display();
 </script>

In my external JavaScript file (my_map.js) has:
    import 'ol/ol.css';
    import Feature from 'ol/Feature';
    import Geolocation from 'ol/Geolocation';
    import Map from 'ol/Map';
    import Point from 'ol/geom/Point';
    import View from 'ol/View';
    import {Circle as CircleStyle, Fill, Stroke, Style} from 'ol/style';
    import {OSM, Vector as VectorSource} from 'ol/source';
    import {Tile as TileLayer, Vector as VectorLayer} from 'ol/layer';
    import VectorContext from 'ol/render/VectorContext';
    import GeoJSON from 'ol/format/GeoJSON';
var my_map = {                       // <-- add this line to declare the object
    display: function () {           // <-- add this line to declare a method

      var view = new View({
        center: ol.proj.fromLonLat([-49.450270, -15.496072]),
        zoom: 6 ,
      });

      var map = new Map({
        layers: [
          new TileLayer({
            source: new OSM(),
          }) ],
        target: 'map',
        view: view,
      });

      new VectorLayer({
        map: map,
        source: new VectorSource()
      });

document.getElementById('produtoid').onchange = function(){
var produtoid = document.getElementById('produtoid').value;
var _token = $('input[name="_token"]').val();
        $.ajax({
                        url:getFeatur,
                        method:"GET",
                        data:{produtoid:produtoid, _token:_token},
                        success:function(result)
                        {

                          var newSource = new ol.source.Vector({
                                        format: new ol.format.GeoJSON(),
                                        features: new GeoJSON().readFeatures(result.dados[0].feature,
                                            {
                                                dataProjection: 'EPSG:26918',
                                                featureProjection: 'EPSG:3857'
                                             }),
                                        });

                          var nwLayer = new ol.layer.Vector({
                              title: 'My Title',
                              source: newSource,
                              visible: true,
                          });
                          map.addLayer(nwLayer); --> I think that here my map has to refresh and show the features.

                      }
                  });
                }

}                                // <-- close the method
};                                   // <-- close the object
export default my_map;               // <-- and export the object


texts: How do you create a loop that assesses if the number is greater than 5. The list [7,2,9,3,1,10,8,6,4,5]?

texts:
x = [7,2,9,3,1,10,8,6,4,5]
for number in x:
      if x > 5:
       print x

This is all I got. New to python.  The output just gives me the list again


texts: Using QGIS, how to save an attribute table as sqlite?

texts:
With 2.14.1 under Windows 10, I am trying to save a .dbf attribute-only table as a .sqlite table.  However, the response is an error.

My procedure: right-click the .dbf in the QGIS table of contents and Save As...  In the subsequent Save Vector Layer As... panel I click the Geometry Type as "No Geometry", as seen in this screenshot:

Here's the resulting error statement.  What am I doing wrong?


texts: Created shapefile layer from selected feature, nothing displaying on QGIS

texts: I selected a feature (watershed) from layer, created a shapefile layer from it, I bring it into QGIS and nothing displays as if there is no data attached to it? Any ideas?


texts: How to install ArcGIS Code Assist in Aptana 3?

texts:
I have downloaded Aptana and ArcGIS Code Assist but I don't know how to "DRAG" the vsDOC into Aptana.

  Aptana 3 Installation Instructions
  Download the zip file that contains the VSDoc file for the version of the ArcGIS API for JavaScript you are using.
  Drag the VSDoc file into your Aptana project.

I drag it into my project workspace and it only shows the .js file content. How do I actually use it?


texts: Moving multiple features on Esri web map at a time?

texts:
Is there is any way to make multi-feature movement on Web AppBuilder map?
I am new to Web AppBuilder and I need to do a move on multiple features at a time.


texts: Google Earth Engine - NASA/NEX Climate Data

texts:
I am working with the NASA/NEX-DCP30 downscaled climate data.

Ideally, I would like to set a threshold for two bands - pr (precipitation) and tasmax (maximum temperature) and filter it to display where those are both met in a given year (e.g. tasmax > 300 AND pr < .0001 for 1995) both visually and ultimately export as a raster and do some charts of how many times it occurred in a year.

Here was my best attempt, maps look nice but it doesn't account for both pr AND tasmax in the same layer...

var collection = ee.ImageCollection('NASA/NEX-DCP30')
    .filterDate('1994-01-01', '1995-12-31')
var maxtemp = collection.select(['tasmax']);
var precipitation = collection.select(['pr']);
Map.addLayer((maxtemp),
    {'min': 286, 'max': 310, palette: '000000, 00FF00'}, 'MaxTemp');
Map.addLayer((precipitation),
    {'min': 0, 'max': .000025, palette: '0000FF, FF0000'}, 'Precipitation');


texts: Why QgsRasterCalculator crashes Qgis unexpectedly?

texts:
I am struggling to sum some rasters using QgsRasterCalculator. When the step "calc.processCalculation()" is executed, QGIS crashes with a  runtime error (R6025
- pure virtual function call) or a dump.

I have a list of raster paths generated in a for loop looking like the following one. Please tell me what do I do wrong and how could I make it work...

Here is my code:

from qgis.core import *
from PyQt4.QtCore import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

s1 = "E:\rcalc\OUTPUT11.sdat"
s2 = "E:\rcalc\OUTPUT22.sdat"
s3 = "E:\rcalc\OUTPUT33.sdat"
list = [s1,s2,s3]

entries = []
formula = ""

for i in range(0, len(list)):
    raster = list[i]
    resp = QgsRasterLayer(raster)
    rast1 = QgsRasterCalculatorEntry()
    rast1.raster = resp
    rast1.ref = "rast"+str(i+1)
    rast1.bandNumber = 1
    entries.append(rast1)
    rast1 = None
    resp = None

    if(i == len(list)-1):
        formula = formula + 'rast'+str(i+1)
    else:
        formula = formula + 'rast'+str(i+1)+' + '

resp = QgsRasterLayer(list[0])

calc = QgsRasterCalculator( formula, "E:\output.tif", 'GTiff', resp.extent(), resp.width(), resp.height(), entries )
calc.processCalculation()


texts: QGIS 3.16 - I want to make by 6"x9" map a 24x36, but all of the text/layers are tiny. Is it possible to automatically scale all of the layers?

texts:
I recently got a freelancing gig and am really excited, but I'm having troubles exporting my map in a different size in QGIS. My client wants a map sized 24"x36", and I have created one that's 6"x9". However, when I go to page properties and change it to 24x36, and drag the map content to fit, the labels and the text are minuscule (I assume 4x smaller).
Is there a way to resize all of the layers ala resizing a jpg in word?  I would think with vectors this would be pretty straightforward. I understand there is a way to do units in map units, but it seems like a big timesink.
6x9 map

24x36 map


texts: Are DBeaver/PostgreSQL Datatype = Polygon real PostGIS spatial Objects?

texts:
I am currently designing a database for my PhD-Project. I want to link a PostgreSQL database to QGIS. I just started in PostgreSQL and am using DBeaver 6.2.2.

When creating my first Test-DB I saw, that I have a multitude of datatype to chose from for a new column, including Polygon, Point and Line. Does that mean I can create a column with this datatype and store my GIS-Geometries in there and is that already a fully functional PostGIS datatype, that is to say could I create e. g. polygons in QGIS from this column?


texts: Calculating sum of geometric length of shapefile in ArcGIS Desktop?

texts:
I have a shapefile with lots of line routes.

I projected it and calculated the length using calculate geometry.

How can I get the sum of all these routes in ArcGIS Desktop?


texts: Executing Select Layer By Attribute in Python Add-in Toolbar (ArcMap)

texts:
I'm having trouble with the button class of my code:

class OKButton(object):
"""Implementation for ok.button (Button)"""
def __init__(self):
    self.enabled = True
    self.checked = False
def onClick(self):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    df = mxd.activeDataFrame
    arcpy.SelectLayerByAttribute_management(layer,"NEW_SELECTION",valueS)
    df.zoomToSelectedFeatures()
    arcpy.RefreshActiveView()

The toolbar allows the user to select a data frame, a layer within it, a field within that layer, and a value within that field. In the above code, "layer" and "valueS" are the global variables referring to the chosen layer and value within one of its fields, respectively. However, when I run the code this error is raised:

ExecuteError: ERROR 000358: Invalid expression
Failed to execute (SelectLayerByAttribute).

I know I need to replace the where_clause with an SQL expression based on the value the user selects (currently the global variable "value"). How could I go about this?


texts: Can ArcMap text formatting support mailto in text elements (ready for Export to PDF)?

texts:
I currently include an email address in a text element created using the Rectangle Text tool and its map gets exported to a PDF file.

Is it possible to use ArcMap text formatting (or any other means) to assign a mailto: hyperlink to that email address so that it persists into the PDF file as something that can be clicked on to start an email?


texts: Mirror / rotate multi links around their starting point

texts:
I have a shp of lines; all the lines have only two vertices
I want all these traits to be "mirrored", or overturned to be clear, around their starting point.
The image explains the concept better; in red I drew the new line I would like to get by hand, while in light blue the original link.
The command must be executed on thousands of links; it is therefore not possible to execute it by line
I know there are already answers on this, but nothing that works on multi links


texts:
I'm trying to select the parcels near a street with web adf for vb.net, does anyone know how to do that? I have a graphics resource in the map resources with the selected streets and I need to do the buffer with all the geometry in that resource and get the parcels around.

Can I use a spatialFilter? How can I get the geometry in the Graphics Resource?