texts: Publishing JPEG2000 layer in GeoServer programmatically (API REST)

texts:
I am trying to publish a JPEG2000 layer as a web map service at GeoServer, firstly, I downloaded and configure the GDAL plugin and all are working fine, so, I can publish manually any JPEG2000 in GeoServer, and after, I want to publish a set of jpeg2000 layers programmatically using the curl command for the API REST of GeoServer using this code:
curl -v -u admin:geoserver -XPUT -H "Content-type: image/JP2" --data-binary "file://C:/Users/pc/Desktop/Images/Jpeg2000/image2.JP2" http://localhost:8080/geoserver/rest/workspaces/cite/coveragestores/image_test/external.JP2

The command does not work with a JPEG2000 file knowing that it works very well for GeoTIFF images. What is the best solution for it to work well?


texts: Convert/export large number of attribute's information shapefiles to txt/csv format using python script?

texts:
I'm using the 2.8.1 Wien QGIS version.
I have a large number (more than 5000) of shapefiles that I want to convert to a csv or txt files to extract, to tables, the atributte information. The shapefiles are located in a directory and I want to keep the same name and directory for the new ones.

Does anyone know how to do this in a python script?


texts: ArcGIS 10.1 Find Identical specify a search order

texts:
I have a large polygon feature class which I know contains duplicate features. I have tried using the Find Identical tool, but is just runs for days without producing a result. I have several different fields that I'm checking against:

Shape
COA
Shape_Length
Shape_Area

The COA field has fewer possible values (thousands vs millions) than the others and duplicates only exist when the COA code is shared. Therefore the logical way to search is to only search for matching shapes in places among fields that have the same COA.

Is there a was to tell ArcGIS 10.1 to search in this manner? I suspect is it currently checking each shape against every other shape in the feature class, which is why it is taking so long.


texts: How to manually reclassify a contiguous region of misclassified pixels?

texts: The attached image is the result of an image segmentation algorithm on high resolution imagery.  The smaller clusters of pixels represent trees and brush and the larger cluster of pixels (marked with a red arrow) is a misclassified farm field.  In ArcGIS 10.1 is there a way to select the misclassified pixels and reclassify them to 0 or NoData?


texts: Find the nearest neighbour for points in one layer to points in another, then measure the distance

texts:
USING QGIS 2.18.19 with python console (planning to upgrade to the newest QGIS after this piece of programming)

I am struggling with finding the correct pyQGIS coding for this.

I have a shapefile containing points marking shelteredHousing and another shapefile containing groceryStores (points).

I wish to identify the nearest grocery store to each fox den point and calculate the distance between them.

This is as far as I have managed so far (and it's taken a long time)

import qgis.utils
layers = iface.mapCanvas().layers()
storesData=0
housingData=0
for layer in layers:
    layerName = layer.name()
    if layerName == "grocerystores":
        groceryStoresLayer = layer
    storesData=1
    if layerName == "shelteredhousing":
        shelteredHousingLayer = layer
    housingData=1
if storesData==0:
    DestinationLayer=iface.addVectorLayer("C:\\destinationpathway\\grocerystores.shp", "groceryStoresLayer", "ogr")
if housingData==0:
    DestinationLayer=iface.addVectorLayer("C:\\destinationpathway"\\shelteredhousing.shp, "shelteredHousingLayer", "ogr")
res = shelteredHousingLayer.dataProvider()

I think I need to have in the statement "from math import sqrt" and then perform a square route on the different layers, but I really wanted to see if someone could give a little help or nudge in the right direction.

Effectively I'm pretty sure I want to do the distance matrix on the data, but I don't know how to plug that in.


texts: Smoothing relatively flat DEM for terrain and hydrologic analysis

texts:
I'm trying to do some terrain and hydrologic analysis from a 1/3 arcsecond DEM that I downloaded from the USGS National Map website. The issue I'm having is that the area in question is central Florida, which is relatively flat. The DEM is a Float64 geotiff, but all the elevation appear to be integer values. What is happening is that in some spots there are large 'flat' areas with a sudden 1 unit drop in elevation. In generating things like slope and water accumulation maps, this is causing clear issues:

(Top: DEM. Middle: slope. Bottom: flow accumulation)

With the slope (middle), there are flat areas with sudden spots of steep slope, when in reality, it should be much more uniform. Similar with the flow accumulation, where there are these sudden areas of low accumulation (white lines) running throughout the map.

Clearly, I need to do some sort of smoothing/interpolation on this map to resolve this. Something simple like a 3x3 or 5x5 sliding window seems like it wouldn't work  because it'll just smooth out the immediate area around the elevation change, but still leave large flat areas that don't necessarily make sense. At the same time, I want to avoid smoothing out the finer details in the more complex areas.

I don't have a really proficient understanding of reconditioning gis data in a meaningful way, so I'm wondering what would be an appropriate method (or methods) of smoothing this data? Primarily, I'm looking for a more theoretical answer than an actual implementation, but if there is also a place I can start looking for an implementation, then that's great as well (I'm on Linux, so no ArcGIS; I use gdal, grass, and R, but can look into other options as well).


texts: How to understand output of ArcGIS System Monitor?

texts:
We had some performance issues with our enterprise ArcGIS server application, in this regard I have installed ArcGIS System Monitor v1.0.9. After installing and configuring the site, they are several variables (below). I would like to know how these variables effect the performance and their definitions. Please throw some light on this.

ARCGIS Percent Processor Time
ARCGIS Working Set GB
ARCGIS Active Working Set GB
ARCGIS Instance Count
ARCGIS Active Instance Count
ARCGIS Idle Instance Count
ARCGIS Virtual GB
ARCGIS Private GB
ARCGIS Page File GB
ARCGIS Page Faults/sec
ArcSOC Percent Processor Time
ArcSOC Working Set GB
ArcSOC Active Working Set GB
ArcSOC Active Instance Count


texts: How to use Mapbox vector tiles as a vector source in OL3, so that labelling will not clip at tile boundaries

texts:
I have been able to load and display mapbox vector tiles in OL3 using a vectorTile source and a vectorTile layer.  The problem with this approach is that labels will clip at the tile boundaries and that features that extend across multiple tiles get multiple labels.  I have increased the guttering around the tiles so that the symbols will not clip, but with potentially large labels, adjusting the guttering will not solve the problem.

I read that vector tiles in GeoJSON format can be used with vector sources, when using a tile loading strategy. Here is my code:

var tileGrid: ol.tilegrid.TileGrid = ol.tilegrid.createXYZ({
            maxZoom: 19
        });

        var vectorSource = new ol.source.Vector({
            format: new ol.format.MVT({
                featureClass: <any>ol.Feature
            }),
            url: function (extent, resolution) {
                var tileCoord = tileGrid.getTileCoordForCoordAndResolution(ol.extent.getCenter(extent), resolution);
                var tileUrl = cadastralTileURL +
                    tileCoord[0].toString() + '/' +
                    tileCoord[1].toString() + '/' +
                    (-tileCoord[2]).toString() + '.pbf';
                return tileUrl;
            },
            strategy: ol.loadingstrategy.tile(tileGrid),
            });

        var vectorLayer: ol.layer.Vector = new ol.layer.Vector({
            source: vectorSource,
            style: this.GetCadastralStyleFunction
        });

        var initCenter: [number, number] = ol.proj.fromLonLat([116, -32]);
        var viewOptions: olx.ViewOptions = <olx.ViewOptions>{
            center: initCenter,
            minZoom: 6,
            maxZoom: 22,
            zoom: 16
        };

        var view: ol.View = new ol.View(viewOptions);
        this.map = new ol.Map({
            layers: [vectorLayer],
            target: 'map',
            view: view
        });

Using the debugger in Chrome I can see the tiles loading properly (the network shows each tile coming down successfully), but nothing shows up on the display. I was hoping to see the polygon features being rendered as whole features, instead of a collection of parts in each tile, and thus getting around the issue of clipped labels.

Can anyone point me in the direction of what I am doing wrong?


texts: Converting mulitlinestring to linestring in PostgreSQL with PostGIS?

texts:
I am using postgres with postgis and I have a set of points which I would like to interpolate with a set of lines.

Now my problem is that the set of lines are of type multilinestring.

Is there a way to convert them to linestring?


texts: Adding group in Legend in QGIS Map composer?

texts:
There is a Map composer with created map. I want to "Add group" to Legend - it appears in the Legend Properties, but not in Legend on the map.

I am using QGIS 2.18.8. This problem was also in 2.16 version. Am I missing something?


texts:
I am exporting a large dataset from PostGIS to GeoPackage and then trying to use it in QGIS and ArcMap. The GeoPackage size is 17Gb and contains around 6 different tables. The biggest table contains around 200million linestring features.

I have created the GeoPackage using ogr2ogr 2.1.3 which automatically adds a spatial index and am using QGIS 2.16.3 and ArcMap 10.4.

When I try and open the GeoPackage in QGIS it takes well over two hours before QGIS becomes responsive and starts to render the features at which point I can add a stylesheet or restrict the zoom ranges.

How can I use this GeoPackage more efficiently?

Do I need to VACUUM ANALYZE it like I would in PostGIS?

Or do I need to somehow add the stylesheet to the GeoPackage to restrict showing features at to high zoom levels.

Extra Details:
Windows Server 2012 64Bit 16Gb RAM


texts: Benefit of adding Metadata descriptions to layers using ArcCatalog?

texts:
I've tried adding Metadata to my primary shapefiles before using ArcCatalog, however, I don't know if this really adds anything professionally to my layers. Is Metadata and Layer description in ArcCatalog really necessary? If so for what benefit?

Currently I run ArcGIS 10.2 and I post finished GIS layers to our County Appraisal District website (http://www.brazoriacad.org/gis-downloads.html), should I be adding Metadata, Layer Descriptions, and thumbnails to our downloadable data or is it more trouble than it's worth?


texts:
I'm trying to achieve a visualization similar to Mapbox GL JS (c.f. Mapbox Vector Tiles from Geoserver 2.11 in Mapbox GL JS), i.e. with pitch and bearing. I'm using the latest OL-Cesium version for this (1.27) and based my code on their examples. This is the most minimal code I can offer (mwe.js):

var style_simple = new ol.style.Style({
    fill: new ol.style.Fill({
        color: '#ADD8E6'
    }),
    stroke: new ol.style.Stroke({
        color: '#3887be',
        width: 1
    })
});

function simpleStyle(feature) {
    return style_simple;
}

var layer = 'opengeo:hessen';
var projection_epsg_no = '900913';
var testLayer = new ol.layer.VectorTile({
        style:simpleStyle,
        source: new ol.source.VectorTile({
            tilePixelRatio: 1, // oversampling when > 1
            tileGrid: ol.tilegrid.createXYZ({maxZoom: 19}),
            format: new ol.format.MVT(),
            url: '/geoserver/gwc/service/tms/1.0.0/' + layer +
                    '@EPSG%3A'+projection_epsg_no+'@pbf/{z}/{x}/{-y}.pbf'
        }),
    visible: true
    });

var osm = new ol.layer.Tile({
        source: new ol.source.OSM()
    });
osm.set('baselayer', true);

var ol2d = new ol.Map({
  layers: [osm, testLayer],
  controls: ol.control.defaults({
    attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
      collapsible: false
    })
  }),
  target: 'map',
  view: new ol.View({
    center: ol.proj.transform([25, 20], 'EPSG:4326', 'EPSG:3857'),
    zoom: 3
  })
});

var ol3d = new olcs.OLCesium({map: ol2d});
var scene = ol3d.getCesiumScene();
var terrainProvider = new Cesium.CesiumTerrainProvider({
  url : '//assets.agi.com/stk-terrain/world'
});
scene.terrainProvider = terrainProvider;

Accompanying HTML:

<?xml version='1.0' encoding='UTF-8'?>
<!DOCTYPE HTML>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    <meta name="robots" content="index, all" />
    <title>olcesium vectore tile layer example</title>
    <link rel="stylesheet" href="../ol.css" type="text/css">

    <input type="button" value="Enable/disable" onclick="javascript:ol3d.setEnabled(!ol3d.getEnabled())" />

When I load the page OL and the vector tile layer load fine. However console is throwing this Error:

Uncaught TypeError: e.ce is not a function

When I click on the "Enable/disable" button the following error is thrown:

Uncaught TypeError: Cannot read property 'setEnabled' of undefined
    at HTMLInputElement.onclick (mwe.html:12)

Is OL-Cesium not capable of using Vector Tile Layers or what could be the problem?


texts: Openlayers and Geoserver. Projections and WMSGetFeatureInfo()

texts:
I use openlayers 3 and Geoserver 2.12.
I'm writing a function to get features by mouse click but in the result I get the features far south-west from the clicked pixel.

I have to notice that my Geoserver's shapes and the group of layers are in EPSG:4326. OL loads my map normally but if try to set the EPSG:4326 projection in the view's declaration it doesn't load it. In debug I can see that view's projection is EPSG:4326 but the map is not shown.

As I know OL sets the EPSG:3857 by default and now I can't change it. So most likely because of the inconsistency of the projections it finds wrong features.

Here is the code:

function initMap() {

var wmsSource = new ol.source.TileWMS({
    url: 'http://178.208.82.80:8082/geoserver/CFD/wms',
    params: {
        'FORMAT': 'image/png',
        'service':'WMS&amp;version = 1.1',
        'LAYERS': 'CFD:CFD',
        'TILED': 'true'
    },
    serverType: 'geoserver',
    transition: 0
});

var layers = [
    new ol.layer.Tile({
        source: wmsSource,
        minResolution: 0
    })
];

var view = new ol.View({
   // projection: 'EPSG:4326', - if uncomment the map will not load
    center: ol.proj.transform([36.587147,50.598299], 'EPSG:4326', 'EPSG:3857'),
    zoom:11,
    minZoom:2,
    maxZoom:18
});

var map = new ol.Map({
    controls: ol.control.defaults().extend([
        new ol.control.ScaleLine({
            units: 'degrees'
        })
    ]),
    layers: layers,
    target: 'map',
    view: view
});

function onMapClick2(evt) {

    var pixel = evt.pixel;
    var bbox = evt.map.getView().calculateExtent(evt.map.getSize());

    var url = 'http://178.208.82.80:8082/geoserver/CFD/wms?' +
        'SERVICE=WMS&' +
        'VERSION=1.3.0&' +
        'REQUEST=GetFeatureInfo&' +
        'FORMAT=text/xml&' +
        'QUERY_LAYERS=CFD%3ACFD&' +
        'version%20%3D%201.1&' +
        'LAYERS=CFD:CFD&INFO_FORMAT=text/xml&' +
        'FEATURE_COUNT=50&'+
        'I=' + parseInt(pixel[1]) +'&'+
        'J=' + parseInt(pixel[0]) +'&'+
        'WIDTH='+ evt.map.getSize()[0] +'&' +
        'HEIGHT='+ evt.map.getSize()[1] + '&' +
        'CRS=EPSG:3857&' +
        'STYLES=&' +
        'BBOX=' + bbox[0] +',' + bbox[1] +',' + bbox[2] +',' + bbox[3];
    fetch(url).then(function(response) {
        return response.text();
    }).then(function(response) {

        var allFeatures = new ol.format.WMSGetFeatureInfo().readFeatures(response);
        document.getElementById('all').innerText = allFeatures.length.toString();
        console.log(url);
        var res = '';
        for (var i = 0; i < allFeatures.length; i++) {
            res = res + ' name = ' + allFeatures[i].getProperties().name +  ' fclass = ' +
                allFeatures[i].getProperties().fclass + ' || ';
        }
        console.log(res);
    });
}

   map.on('singleclick',onMapClick2);
}


texts: Turn bathymetry into elevation in arcGIS?

texts: My raster is in depth of bathymetry however when this is shown in 3D it reverses the elevation making the deepest points the highest in elevation and the shallowest points the deepest. This makes a trough look like a mountain when displayed in ArcScene. The data is in distance from the ship not in elevation from the sea bed, is it possible to reverse the data? How would I solve this?


texts: How to create efficient database for routing?

texts:
I created a database of many important features and locations of my city, and I want to find shortest/fastest path between two locations. Now, I'm calculating the route by using the A* algorithm in Java for Android device, but it takes too much time to calculate and draw a simple route.

Is there any efficient way to store and retrieve the road data from database for route calculation?


texts: Seeking dataset locating points of interest such as fire stations in Canada?

texts:
Is there a dataset locating points of interest such as fire stations in Canada?

I tried using the one at http://www.mapcruzin.com/ but the accuracy was not correct.


texts: Can QGIS import Multipatch geometry type?

texts: can QGIS import multipatch 3d geometry type? What is the data structure of multipatch how computer record and draws it ? Is it solid or not??


texts: I have been doing some research and found out that i can change the attributions of TileWMS, OSM but i have found little about BingMaps. Is it possible to change the attributions of BingMaps? and how do you go about it?


texts: Surface attribute only calculated for actual points and not interpolated

texts: When I run the surface hillshade, slope and aspect in the Arc Toolbox using my raster surface it appears to be calculating the attribute only for the data points and not for the interpolated areas between. I loaded the points as a XYZ file then used Feature to Raster to create the surface.


texts: How to split summary address received as response of ordnance survey address base

texts:
I am using an API to fetch search response from ordinance survey address base, the result includes a summary address field which contains full address in one line.

The result can be ambiguous as summary address can contain any type of details in it. An address can have 4 or 5 comma separated fields or sometimes even 8 or 9 comma separated fields.
example 1 : "results": [{
    "address": [{
      "uprn": xxxxxxxx,
      "summary_address": "2 Wakefield, Edeinnburh, EH9 6TL",
      "matchscore": 100
    }

example 2: "results": [{
    "address": [{
      "uprn": xxxxxxxxx,
      "summary_address": "Domino Piza, Unit 10, 100 sample Road, Bishopbrig, G44 2AB",
      "matchscore": 100
    }

I have no clue how to split these fields properly and identify them. Is there any fixed format for this?
I want to split address and identify each field either as house number, house name, road, locality, town etc

Note: the address used in above examples are fake and do not correspond to any real data


texts: Rotation of a cropped portion of a raster map

texts:
We create printed maps from a GeoTIFF file projected in LCC, which is a common projection for aviation maps. Part of our process includes customization of the raster image by selecting a region and adding other design elements. The final part of the process allows for the customer to select desired airports on the map via a custom leaflet application showing all the airports via a GeoJSON file.
We would like to have the tiles map image the customer sees in the Leaflet application be similar to the way the printed map will appear to the customer. This means they will stay in LCC and we crop and raster tile the actual artwork for use in leaflet. We have this all worked out just fine with tiling the resulting image and using leaflet and proj4leaflet. The issue we have is sometimes would like to rotate the crop so that the raster image does not end up aligning with the LCC projection.
Shown here is a copy of my main image where you can see the inside image which will make up our final printed artwork. The inlayed image aligns with the original image, but the cropped area is rotated about 13 degrees.

I would like to output inset cropped image, rotated so that it appears as it would appear hanging on the wall. This would involve rotating the image about 13 degrees.

If it were not rotated, we easily calculate a new TFW and include the northing and easting components by calculating the pixel size and the amount of the crop. This all works out perfectly.
However, with it being rotated, we need to add the rotation components to the TFW file and I have been unsuccessful in trying to figure that out.
The end results should be a map which appears to be a cropped version of the main map, projected in LCC, but rotated 13 degrees so that it appears online as it would appear hanging on the wall.
Things I've Tried
I've studied the page here: https://en.wikipedia.org/wiki/World_file and worked out that I need to do some trig on the pixel values.
The other completely unknown is if Leaflet and proj4leaflet can even handle these translations? Answer -- YES, I think I can use the Leaflet.Transformation function for this once I figure out the transformation to use for the rotation.
Update
Here is a concrete example of what I am trying to do along with the files.
The projection of these are
'+proj=lcc +lat_0=38 +lon_0=-100 +lat_1=45 +lat_2=33 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
Here is a portion of the the map I'm working with along with the accompanying pgw file.

Here is that same image rotated 45 degrees (around the center of the image)

What should the corresponding world file be for the rotated 45 image and how do you calculate each field in that file.


texts: Standard user interface for QGIS plugins

texts:
QGIS plugins found in "processing" all have the same user interface.
In particular, to open a file, the selector either show the list of already downloaded files, or have a "..." button to browse and load a file, which will replace the list of files.

I could (painfully) rewrite this code, but I'm wondering if there are some ready-made code to generate this kind of "standard" user interface?


texts: Configuring QGIS to show only toolbars wanted?

texts:
I friend of mine asked me a question the other day that I didn't knew the answer to. In his organization, they are using QGIS for non GIS-experienced people. He would like to configure QGIS to only show the tools/toolbar that those non-experienced people need to see. A common complaint is that QGIS looks too so difficult with all toolbars. QGIS is installed to this users on a server and he would like to control this so that people don't need to change the toolbars on every computer that are using QGIS.

How can I solve this kind of problem?


texts:
I am trying to align two shapefiles (LA County freeways and LA County city boundaries) to see which freeways run in each city.

In QGIS, the shapefiles don't align. The freeway layer has a warning question mark saying "no CRS set for layer," however once I set the CRS to match the city boundaries layers (WGS 84 epsg:4326), nothing changes. I "save feature as" shapefile with the same CRS and the imported result is still misaligned with the boundaries.

I tried to export the file again, setting the extent to match the city boundaries, but then it just exports with 0 features.

How can I get the files to align? Specifically, how do I get a file with no CRS to align with another if simply setting the CRS and re-importing does not work?


texts: Changing default parser in ArcGIS Field Calculator to Python?

texts:
Is there a way to change the default parser in the Field Calculator to Python?

So annoying to always change now that I'm trying to move away from using VBScript...

I'm using ArcGIS 10.1 for Desktop.


texts: ArcGIS Server JSON to GeoJSON does not always work

texts:
I'm using GDAL 1.10.1 for reading ArcGIS Server responses with automatic conversion to GeoJSON format as described in Converting ArcGIS Server JSON to GeoJSON? post.

However, it looks like some of the Feature Servers published by ArcGIS return empty results when queried using the above technique:
http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/BloomfieldHillsMichigan/LandusePlanning/FeatureServer/2/query?where=objectid+%3D+objectid&outfields=&f=pjson

Could someone explain why the technique is not working for some data?


texts:
I want to export my classified image to drive in GEE, for presenting in a document.
I have used this code.

Export.image.toDrive({
  image: classified,
  description: 'classified_basin_with_points',
  scale: 30,
  region: geometry,
  fileFormat: 'GeoTIFF',
});

But the exported image showed no color whereas my cassified layer have 4 clssifying color on it. How can i get a colorful image that will cover my classifying layer color?


texts: Where can I find GDAL Python bindings for Python 3.4?

texts: I am using python 3.4 for development and I need GDAL python bindings for python 3.4. I have already gone through this site, but I can only find up to python 3.3?


texts:
I'm developing a simple app with ESRI iOS SDK. I set up esri-sdk for iOS with CocoaPods, but CocoaPods installs all frameworks in project. In Xcode theestimated project size is about 250 MB.

I'm curious when I run ipa file, how many megabytes the ipa file will be on the device and then how can I decrease app size?


texts: Field length in geopandas shapefile output

texts:
I am converting an excel sheet to .shp via pandas i.e. I have a dataframe with 150+ columns and ~900 records (increasing by about 50 every day), that I convert to geodataframe and export to shp:
output_data.to_file('myfilename.shp', driver='ESRI Shapefile')
Saving this data to shp in Global Mapper creates a 3.5MB dbf
In QGIS, a 14MB dbf
Using Geopandas 6.5MB dbf
I've seen that the geodataframe has 15 decimal places for every numeric field, including whole numbers. Same thing for QGIS. Global Mapper's field lengths are equal to the lengthiest field in each column.
How can I keep geopandas from adding 15 decimals? I assume this is what bloats the filesize since I refactored some fields in QGIS and got a much smaller output dbf. The fields are also shorter in pandas df, before converting to geodf.
Edit:
Couldn't find a way to use GDAL RESIZE before saving as shp.. What I ended up doing:
export to shp using geopandas
output_data.to_file('myfilename.shp', driver='ESRI Shapefile')
resize using gdal/ogr
dbf = ogr.Open('myfilename.shp', 1)
dbf.ExecuteSQL('RESIZE myfilename', dialect='sqlite')
This took a 47MB dbf down to 27MB, however only shortened string fields.


texts: Use Mapbox featureLayer to create Marker Clusters

texts:
I'm trying to use this tutorial (https://www.mapbox.com/mapbox.js/example/v1.0.0/markercluster-with-mapbox-data/) to create marker clusters to style a Mapbox featureLayer containing my data. This code works to display my data:

var featureLayer = L.mapbox.featureLayer(rodents1,{
        pointToLayer: function (feature,latlng){
            return L.circleMarker(latlng,markerStyle);
        },
    }).addTo(map);

But when I try to rework my code according to the instructions, my data goes missing from my map. Here is the code according to the instructions:

// Since featureLayer is an asynchronous method, we use the `.on('ready'`
// call to only use its marker data once we know it is actually loaded.
L.mapbox.featureLayer('rodents1').on('ready', function(e) {
    // The clusterGroup gets each marker in the group added to it
    // once loaded, and then is added to the map
    var clusterGroup = new L.MarkerClusterGroup();
    e.target.eachLayer(function(layer) {
        clusterGroup.addLayer(layer);
    });
    map.addLayer(clusterGroup);

I've included links to the Leaflet/markercluster files in the head of my html file.

Am I missing something?


texts:
I have a geotiff (RGB) 8-bit per channel file that was classsified to four classes. Theoreitcally I could save the file as a 3-bit and make the file significantly smaller. For my application files size is very important. Yet from what I can tell from the GDAL documentation http://www.gdal.org/frmt_gtiff.html I think am out of luck. Since each pixel is a stand-alone, aggregating the pixels and converting them into a polygon does not seem it would do the trick. For my processes I have been using GDAL API with Python 3.5.
Does anyone have any thoughts on how to take advntage of the classified data to convert to s smaller file size?


texts:
I am working with ArcGIS Network Analyst (10.5.1).  I created a network in SDE and everything worked fine.  The network is versioned on SDE.

I updated the source line features of the network in a geodatabase on my local drive by adding new line segments, rebuilt, tested it, and it still performs as expected.

I then used Delete Rows and Append tools to load the updated line features into the network on SDE and rebuilt it.  Routing and other analysis still work, but now the performance is much slower.  A routing analysis that used to take a couple seconds to solve now takes close to a minute.

Has anyone else experienced this? Anyone have any thoughts what can be done to optimize or speed up a network stored in SDE?

I suspect if I delete and recreate the network in SDE it would work OK again, but I would like to avoid going through that whole process every time I need to update it


texts: Converting geometry to geography using org.geotools.geometry.jts.JTS?

texts:
I'm trying to convert com.vividsolutions.jts.geom.Geometry to geography, using the GeoTools package - org.geotools.geometry.jts.JTS. This is my code:

public Geometry getGeography() {
    Geometry geo = null;
    try {
        JTS.toGeographic(this.geo, CRS.decode("EPSG:" + Integer.toString(this.geo.getSRID()), true));
    } catch (TransformException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    } catch (NoSuchAuthorityCodeException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    } catch (FactoryException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
    return geo;
}

I'm getting:

org.opengis.referencing.NoSuchAuthorityCodeException: No code "EPSG:0" from authority "EPSG" found for object of type "EngineeringCRS".

Can anybody assist on why the this.geo.getSRID() return 0?


texts: Roles for uploading layers in geoexplorer

texts:
What kind of users I have to create for uploading layers? I tried with users that are not admins and they have specific role with write access on workspace, but it works only once or twice. What I need is following: many "public" users with anonymous login and access to public workspace (stores and layers) and several power users which can see all in public workspace but they have their common private workspace for writing. Private workspace is for uploading their shapefiles and geotifs, and that workspace first users don't see. So, is it possible to have non-admin users for that feature in geoexplorer?
OS is Windows 2008 Server, I don't know is it necessary to give some permissions on tmp\upload and private folders?


texts: Buffering around feature except where barrier prevents using QGIS

texts:
Is is possible to grow a buffer around a feature, except where another feature stops the expansion using QGIS?

There is an existing similar question for ArcGIS Desktop: Buffering with physical barrier using ArcGIS Desktop?

My concrete task is finding cycleways suitable for children, where one criterium is "not dangerously close to major road". I can do this with a buffer around major roads, and removing cycleways within this buffer. However, if some suitable barrier exists, I would like the buffer/"danger zone" to stop there.

Simply making another buffer around the barrier and subtracting from the "danger zone" buffer around the road won't do, as that could also erroneously remove cycleways on the "dangerous" size of the barrier from the "danger zone". I was thinking of somehow finding which is the "safe side" of a barrier, and make a "safe zone" feature from that to cut from the "danger zone" buffer.


texts:
I'm trying to migrate my desktop GIS application from Bentley Map to QGIS. Everything seams to be transitioning great, expect the rotated points. The rotation angle they were saved in Bentley was in radians. I have tried using the mathematical formula (degrees = radians * 180 / pi) to convert them, but It doesn't get me the right angle.

Is there another way to do It?


texts: Trouble display Chart in PDF exported from ArcGIS

texts:
I am using ArcGIS 10.1 and we have MS Office 2013.  I am not sure if this has a lot to do with the MS Office 2013 that ArcGIS 10.1 doesn't read...So I am sharing this with you all...

When I set it up a chart in MS office 2013. and then copy and paste to ArcGIS 10.1 it will show the chart there but when you export to any fomrat  the chart will not display.

Has anyone had any issue with it ?  Is there a way around to fix this problem. The only way is I have to do is save it in Power Point and then save it as a TIFF file and then inset it in ArcGIS .  It is okay...

The first one you see is in ArcGIS Layout View. The second one you see is export in PDF .


texts: Using Ripley's K Function in ArcGIS for coral density along a transect?

texts:
I have some data from transects taken in a mangrove system. I want to analyze the spatial distribution /some sort of spatial statistics of how clumped/dispersed corals are along this linear transect. We have GPS coordinates for each coral, as well as the measured distances between each coral.
Does anyone have knowledge of how to conduct a spatial analysis of points along one line, using possibly Ripley's K or another function?


texts: How to install Sextante in QGIS 1.8 without Plugin Installer?

texts:
How can I install sextante (and where can I download it) if the Plugin Installer doesn't works?

As I work for a corporation, I have some problems conecting to the official repository, and I cannot solve them ...


texts: What is the difference between Vector and Raster data models?

texts:
From: http://support.esri.com/en/knowledgebase/GISDictionary/term/vector%20data%20model

vector data model: [data models] A representation of the world using points, lines, and polygons. Vector models are useful for storing data that has discrete boundaries, such as country borders, land parcels, and streets.
raster data model: [data models] A representation of the world as a surface divided into a regular grid of cells. Raster models are useful for storing data that varies continuously, as in an aerial photograph, a satellite image, a surface of chemical concentrations, or an elevation surface.

All I have understood from the above is that both vector and raster data constitute of "latitudes and longitudes", only. The difference is in the way they are displayed.
Latitudes and Longitudes in Vector data are displayed in the form of lines, points, etc.
Latitudes and Longitudes in Raster data are displayed in the form of closed shapes where each pixel has a particular latitude and longitude associated with it.
Is my understanding correct?


texts: Cloning multiple layers within group using PyQGIS

texts:
I have a problem cloning more than one layer within a group. I am not proficient in PyQGIS by any stretch of the imagination. My task is to effectively duplicate all raster layers which names end with "_DEM", and append " copy" to the newly formed duplicates or clones. My code snippet for doing only one cloning is as follows:
names = [layer.name() for layer in QgsProject.instance().mapLayers().values()]
substring = "_DEM"
strings_with_substring = [name for name in names if substring in name]
print(strings_with_substring)

DEM = QgsProject.instance().mapLayersByName(strings_with_substring[0])[0]
duplicated = DEM.clone()
QgsProject.instance().addMapLayer(duplicated, False)
group.insertChildNode(3, QgsLayerTreeLayer(duplicated))

dname = duplicated.name()
duplicated.setName(dname + " copy")

Is there a way to construct a function or a loop to iterate this piece of code through every instance of the layer list i.e. strings_with_substring?


texts: Count number of monthly water pixels from JRC Water History Collection using Google Earth Engine - Rgee

texts:
I'm trying count the number of monthly water pixels from JRC Water History Collection using Google Earth Engine - Rgee, over a set of regions.
https://developers.google.com/earth-engine/datasets/catalog/JRC_GSW1_2_MonthlyHistory
This is what i am trying to achieve:

Select the image collection and filter by dates for the time series.

Select water only pixels and create masks containing only water pixels for each monthly image. Note this dataset uses a Bitmask:
Bits 0-1: Water detection
0: No data
1: Not water
2: Water

Count the number of water pixels in each region for each monthly image (create regional time-series)

I think something is going wrong in step 2 - as some pixel counts are missing, i'm not confident the results are correct, can anyone help?
library(rgee)
library(mapview)
library(sf)

ee_Initialize

wnf_shapes <- st_read("./wnf_shapes.shp")

# Step 1
filteredJRCDataset <- ee$ImageCollection("JRC/GSW1_2/MonthlyHistory")$
  filterDate("2006-01-01", "2006-12-31")

# Step 2
surface_water <- filteredJRCDataset$map(function (image) {
  waterImage <-image$select('water')$eq(2)
  image = image$updateMask(waterImage)

  })

# Step 3
ee_sw <- ee_extract(x = surface_water, y = wnf_shapes,  scale = 100, fun = ee$Reducer$count(), sf = FALSE)

colnames(ee_sw) <- sprintf("%02d", 1:12)
ee_sw$id <- wnf_shapes$id

link to shapes files - https://drive.google.com/file/d/1oWJ_ZpEQ4bEYr7R73oOXrQc9UhOH_oCB/view?usp=sharing
I could also do this using Javascript or  Python APIs for anyone not using Rgee.


texts: qgis-processing: possible to specify variable number of output layers?

texts:
Qgis-processing advises that output layers are specified with syntax such as

##foo=output vector

The manual advises against manually loading output layers (in your script code) as QGIS will do this automatically based on the above code.

What if the number of output layers is not known until runtime?  Is there a way to do this besides loading them manually and specifying ##nomodeler?


texts: How to create a polygon by reading the geometry of another polygon layer using ArcPy?

texts:
Basically, I am trying to generate polygon geometries based on the selection of another polygon layer. These polygons will be copied into a new empty schema. I avoiding using the Append and Copy operations since they won't meet all the requirements needed for my output. I also found the cursors help me achieving the desired result way faster. As a result, I need to print the polygon coordinates as one single line text string instead of multi-line. This output will be used as a variable to generate the polygon later in the code.

for row in arcpy.da.SearchCursor("countAddrPnts_lyr", ["OID@", "SHAPE@", "SITUS_STREET_NUMBER_END2", "SITUS_STREET_NUMBER2"]):
x = 0
for part in row[1]:
    for pnt in part:
        if pnt:
            pntArray = ("arcpy.Point(" + "{}, {}".format(pnt.X, pnt.Y) + "),")
            print pntArray
    x += 1

The current output looks like this:

arcpy.Point(907781.079669, 628629.095669),
arcpy.Point(907782.765669, 628579.116669),
arcpy.Point(907649.596669, 628574.910669),
arcpy.Point(907647.068669, 628649.878669),
arcpy.Point(907780.236669, 628654.092669),
arcpy.Point(907781.079669, 628629.095669),

I need it to look as follows:

arcpy.Point(907781.079669, 628629.095669), arcpy.Point(907782.765669, 628579.116669), arcpy.Point(907649.596669, 628574.910669), arcpy.Point(907647.068669, 628649.878669), arcpy.Point(907780.236669, 628654.092669), arcpy.Point(907781.079669, 628629.095669)

Maybe my approach is not the correct one, so any suggestions will be welcome.

BTW, I am not a programmer per se therefore, my code will not be the prettiest but it is doing the job.

In this other post you can download the data in GDB format, as well as the py script.


texts: Create PDF on server side by passing layer from flex client

texts:
I want to pass the user selected layers from client

User can  select and change the order of layers as populated from different map services hosted on my arcgis server
User can set the transparency for each layer

I want to create the web services for that I found some examples but i due to some reasons I don't want to use SOE


texts:
I'd like to know how to select duplicates records except one, with the "select by attribute".
I've attached a picture to explain better the output needed.


texts: Extracting lateral averages using ArcGIS Desktop?

texts:
I want to do the following: to calculate the lateral average value of the surface aspect of a glacier (i.e. the average along a specific contour, the greenish lines).

Can someone help me or say how this can be done?


texts: Is there any random data generator supporting geometry in PostGIS?

texts:
I found few random data generators for PostgreSQL (such as this one), but none supporting PostGIS geometry type, other than point. Making points would be quite easy, but for lines with multiple nodes and for polygons (I need more complex polygons than simple rectangles) it would be hell to write the function myself. I would prefer using others' tools (free and opensource, if possible) instead of making them myself.

EDIT: there are algorithms generating random points within a polygon. For polygons, I would prefer a different approach: to generate quasi-centroids first and then make polygons around them. I have some idea how to do it (random starting direction or random angle from previous node + random distance from quasi-centroid). Although if I found a tool generating polygons from scratch, making a centroid or a random point within would be easy. If no one points to any such a tool in a few days, I'll make it myself and post the code as an answer.

I don't need complex line generator so much now, so I'll either make it later, or wait for some other answers. Working complex line generators will get at least +1 from me even if I found I won't really need them, and if combined by some good complex polygon generation they should merit acceptation.


texts: What is the best geocoding Plugin to use in QGIS?

texts:
I have to transfer around 1200 points into my GIS project. Therefore I have given a .csv-file with city names (partly missspelled or unidentifiable).
All of the locations are in Germany (as I read that some geocoding Plugins are for specific places only).
What I need the Plugin to do is read out the .csv and create a georeference point layer in each of the (readable) cities.


texts:
I have Desktop 10.1 that is accessing data from sde geodatabase.
I select some features on the map, and try to do "Create layer from selected feature"
but nothing happen, and no error message.
the operation was successful when I had sde geodatabase 9.3, but since I migrated sde geodatabase to 10.1 and I cannot do that.

The RDBMS is oracle.

I check for errors in
c:\users{my user}\AppData\local\Esri\
and nothing related to any error

or this location
c:\users{my user}\AppData\roaming\Esri\

Where are the selection layers saved?
are they saved in the map document?


texts:
I have two shapefiles showing LULC of the same area for two different years; 1990 and 2015. Number of Polygons in two files are different but the number of class in both are same.

Is there any way to do the change detection study using these vector files?

I have ArcMap 10 and eCognition to use? Suggestion of any other tool is welcome.


texts:
How can I automatically center the map view on the drawn GeoJSON?

My GeoJSON can either contain a polygon or multiPolygon feature.
I'm drawing the following GeoJSON onto my Leaflet map:

let parsedJson = JSON.parse(myString.geoJson);

let geoJSON = L.geoJson(parsedJson, {
     ...
}).addTo(myMap)

I tried fitBounds on the map object, but both center the map out of bounds. Zoom: I'd like to set the view to the bounds of my polygons.

Edit: I realized, the coordinates returned by getCenter() are for example LatLng(175.836357, -36.649876). When I switch them to LatLng(-36.649876, 175.836357) the center point is correctly drawn.

If I'm not mistaken, LatLong ranges are:

             North (+90)
               |
(-180) West ———+——— East (+180)
               |
             South (-90)

So the point returned by Leaflet should be correct and mine should not make any sense - I don't understand why it's correctly drawn when flipped. Could be the reason why the polygon.getBounds() LatLong's are also flipped.

Workaround: I found this working when switching the latLng coordinates of the bounds calculated by Leaflet. Sadly this turn the topics question to "Why are my coordinates flipped?" and leaves me with the workaround of just flipping them;

let coordinates = parsedJson.geometry.coordinates;
let polygon = L.polygon(coordinates);
let bounds = polygon.getBounds();
let southWest = bounds.getSouthWest();
let northEast = bounds.getNorthEast();
let cSouthWest = L.latLng(southWest.lng, southWest.lat);
let cNortEast = L.latLng(northEast.lng, northEast.lat);
let newBounds = L.latLngBounds(cSouthWest, cNortEast);

myMap.fitBounds(newBounds);


texts:
I am using ArcMap 10.5.1

I have been using the Clip tool to cut out the amount of a given habitat type found within individual buffers around individual samples. After successfully executing this command many times and receiving outputs, the Clip tool suddenly started to produce an error message:

  Error000210: Cannot create output

despite my having made zero changes between successful/unsuccessful executions of the command.

I then switched to the Intersect tool which worked fine for a while, until it too started to produce a warning message:

  Warning000117: Warning empty output generated

again without my having made any changes.

I have tried to "Check/Repair Geometry", as well as recreating/reloading my shapefiles, but to no avail. Also, the errors now occur on all samples, including the ones that originally worked.

Any suggestions?


texts: (Officially retired) ESRI Silverlight API

texts:
I have a client with a web GIS app based on the now retired ESRI Silverlight API.

Eventually we will be replacing this app with a javascript based alternative, but for now we want to be able to load in layers from CSV, KML and WMS which apparently is possible via an additional assembly (ESRI.ArcGIS.Client.Toolkit.DataSources)
i.e.

https://developers.arcgis.com/silverlight/api-reference/#ESRI.ArcGIS.Client.Toolkit.DataSources~ESRI.ArcGIS.Client.Toolkit.DataSources.CsvLayer.html

This is obviously not available officially via ESRI anymore, and the client cannot find the original download of the API. All we have is a .NET app with the relevant assemblies (but not the one I need) added as a reference.

Does anyone know where these retired API downloads are still available?


texts: From DWG to Shapefile: How to extract the channels length of water courses?

texts:
I have a DWG file with the drainage system of a small town. I need to open it on ArcGIS and get the total length of the channels that are part of one of the watersheds in the city.

When I convert the DWG to shapefiles, the attribute table on ArcGIS does not show the length of each channel?

I am doing anything wrong? What is the correct way to convert a DWG file to shapefile and be able to also export the geometric features of the watershed?


texts: Calculate Field - getting the code value but not the description

texts:
If STREETNM value are blank or _ then I want to replace them with the value of ST_DESIG. My code below returned the code value of ST_DESIG but not the description.

How can I get the description value of the ST_DESIG. For example, "Divided Major Highway - Class II" NOT 1.

import arcpy
Streets_Centerlines = r"G:\xStreetNew\Test.gdb\Streets_Centerlines"
expression = "replace(!STREETNM!)"
codeblock = """
def replace(STREETNM,ST_DESIG):
 if STREETNM.lstrip().rstrip() == '' or STREETNM.lstrip().rstrip() ==    '_':
      return ST_DESIG
 else:
      return STREETNM
"""
arcpy.CalculateField_management(Streets_Centerlines, "STREETNM", "replace(!STREETNM!,!ST_DESIG!)", "PYTHON", codeblock)


texts:
When I process a raster in QGIS it creates a second band. It appears that the band has no useful data. I like to use the Raster Colours Tool and it will not accept a raster with anything more than one band.
How do I remove the extra band?


texts:
I have a .nc file  that you can download here. This is satellite imagery.

I want to re-project this and save it to a new .nc file, without saving to raster. Just NetCDF to NetCDF only. Is this possible?

I am getting an error saying "ERROR 4: `test.nc' not recognized as a supported file format."

The code I have tried is :

gdalwarp -t_srs EPSG:3857 test.nc test-projected.nc

I am not a professional with this, but should I even be using gdalwarp to do this? (maybe gdaltransform?) I only want to change the projection - nothing else at this time.


texts: Joining two polygons that have information in same spatial location using ArcMap

texts:
I am using ArcMap 10.6.1.

I have two shapefiles, one that has human disturbance and one that is fishnet grid and every grid has a unique value. I used the fishnet tool on my human disturbance shapefile, to create grids of the same size so that for each grid I can say how much area within a grid is human disturbance.

How do I join these two shapefiles, so that I can have the unique value from one shapefile and the human disturbance of the other all in the same polygon?

I have tried spatial join and merge, and they have both returned null values.


texts: Accuracy assessment using polygons as reference data as against point data?

texts:
Is it possible to do accuracy assessment using polygons as reference data as against point data?

I have field-identified polygons as data for training and reference (from a forest-agri mosaic) that have been identified from field. The regular post-classification accuracy assessment method is to use random points and check background image. But instead can I use my field polygons as reference data, instead of points. Problem being, I cannot correctly interpret the land use on the background image or even on googleearth since it is complex forest-agri mosaic. All I have is my reference data which are polygons.

I could convert each polygon into a point and do it, but that defeats the purpose, I feel. the idea is to try to get the software to recognize all pixels in every reference polygon as a single class. Can this be done? If yes, how? I am using ERDAS 2014, but I also notice the QGIS 2.6 plugin for accuracy assessment.


texts: QGIS 2.2 saving layer as CSV with GEOMETRY=AS_WKT not working

texts:
Just upgraded to QGIS 2.2 and it looks like saving a layer as CSV doesn't support GEOMETRY=AS_WKT anymore. The new Layer Options portion of the save dialog is cool but it looks like GEOMETRY and LINEFORMAT for CSV have some swapped options (CRLF and AS_WKT are in opposite option lists). I tried custom layer options:

GEOMETRY=AS_WKT
SEPARATOR=SEMICOLON
LINEFORMAT=CRLF

Those worked under QGIS 2.0 but don't appear to work with QGIS 2.2.


texts: Rasterize does not rasterize field with the value 0

texts:
I have a shapefile layer that contains the field fraction which varies from 0 to 0.05.
DN       sum       count       fraction
1        0         90          0
2        5         1070        0.0047
3        2         39          0.0051
4        1         164         0.0061
...      ...       ...         ...
1615     0         90          0

I want to rasterize the fraction field with the GDAL function Rasterize in qgis.
I set the width and height to 349 & 277, respectively. The map extent is set to another loaded map in the canvas. But, when I execute the rasterize function it excludes all the values of 0 in the plot. Below is my work flow!
How can force the function to include the values of 0 in my raster plot?

The output of my 'struggle' is given here. Purple indicates what shape it should have resembled, black indicates my rasterized file.


texts: python 2.4 replica synchronize changes get system error (0)

texts:
I have python scripts to synchronize changes between 2 geodatabases, the script working well before and after I remove 2 feature classes from one replica, and create a new replica using these 2 feature classes, all my replica cannot be synchronized using the following script:

import sys, string, os, arcgisscripting, smtplib

# Create the Geoprocessor object
gp = arcgisscripting.create()

# Set the necessary product code
gp.SetProduct("ArcInfo")

gp.AddToolbox("C:/Program Files (x86)/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")

# Local variables...
replica_sde__2_ = "C:\\Connections\\replica.sde"
user_prod_sde__2_ = "C:\\Connections\\user_prod.sde"

# Process: Synchronize Changes...

gp.SynchronizeChanges_management(replica_sde__2_, "DBO.Replica1", user_prod_sde__2_, "FROM_GEODATABASE2_TO_1", "IN_FAVOR_OF_GDB2", "", "")

I got system error (0) and I failed to find any detail for this error.

I am using python 2.4 and ArcGIS desktop 9.2.

Anybody has the experience on this?


texts: Google Earth Automated Imagery Downloading for different locations

texts: Is there any method for Google Earth to Automate the process where we define cities and areas and we receive printable files from the most interesting spots of the cities?


texts: Joining Two Fields of Differing Cases in QGIS

texts:
I am attempting to conduct table joins using a common ID field in QGIS. The issue it that the ID field in one table contains letters in the upper case and the other contains them in the lower case. For example, table 1 would have the ID SAM and table 2 has the ID sam. For my purposed a differing case does not make two ID unique if the string is the same otherwise.

I have tried converting the IDs the the same case using the UPPER() and LOWER() function in the field calculator, but due to the size of some of the datasets this can be time consuming and rather inefficient. Is there a way to conduct table joins in QGIS in a manner that is not case sensitive?


texts: Select layer by attribute - SQL command pending on two columns

texts:
I have a layer, with two columns called NAME and TYPE.
A polygon in the layer has only one NAME, but can contain multiple TYPE of vegetation.
I want to select different Types of vegetation, where the name is the same.

Example:

I want to select all the polygons that contain Evergreen AND Dipterocarp.

From the pic the polygon named Doi chiang dao should be selected.

I have tried using the "Select Layer by Attribute" in ArcMap, but I cannot figure out how to create an SQL statement that can do the above-mentioned operation.


texts: How to use raw_input() in QGIS python console?

texts:
I'm hoping to use raw_input() inside a loop to prompt the user when s/he wants to continue running the script. However, I'm getting an error:

Traceback (most recent call last):
   File "<input>", line 1, in <module>
EOFError: EOF when reading a line

Is this function supported in qgis' python console? If not, do you have any suggestion?

I don't need to read the user's input. I just need the loop to pause for a while to let the user do some other stuff, then when ready, to resume working again.


texts: Working out percentage of polygon covering another using PostGIS

texts:
I have 2 polygon tables, polytable_a and polytable_b. a has 1000s of features and 100s of unique names (a1, a2, a3, a4.... etc) assigned to the polygons. b also has 250,000 features but only 3 unique names (b1, b2, b3).
No features within table a overlaps other features in table a. The same is the case for table b.
I'd like to work out, for each unique poly name in table a, what % of it's own area is covered by the various shapes within table b?
ie:
polytable_a name    b1%   b2%   b3%
a1                  40    5     0
a2                  0     20    5
a3                  0     0     5
a4                  0     0     25

bm13563's answer below works great for small poly tables, but the query runs endlessly for these 2 large tables I'm using. I need a solution that will work on 2 tables with 100s of thousands of rows.


texts: Distance and Functions in Google Earth Engine: Producing function that counts points at a given distance each year

texts:
I have tile centroids tileCentroids that cover a polygon in Sierra Leone. I also have a feature collection of points of interest in Sierra Leone saved as camps. The points are only present in some years, but not in others. This information is stored as a series of binary variables open_2001 - open_2012 in camps.
I am trying to produce 11*20 variables indicating how many of the points were present in each year (2011-2012) that were 1 km away, 2 km away, ...20 km away from each tile centroid. I have a function that performs this for each individual years. An abbreviated version of this function, which only gives us the number of points in 2003 that are 1 km from the tile centroid, is:
var spatialFilter = ee.Filter.withinDistance({
  distance: 80000000,
  leftField: '.geo',
  rightField: '.geo',
  maxError: 10
});

var joined_3 = ee.Join.saveAll({
  matchesKey: 'neighbors',
  measureKey: 'distance',
  ordering: 'distance'
}).apply({
  primary: tileCentroids,
  secondary: camps.filterMetadata('open_2003', 'equals', 1),
  condition: spatialFilter
});

var withNearestDist= joined_3.map(function(f) {
  var nearestDist = ee.Feature(ee.List(f.get('neighbors')).get(0))
      .get('distance');
  return f.set('nearestDist', nearestDist);
});

var pointdata_withbandcounts_3 = withNearestDist.map(function(pt) {
  var neighs = ee.List(pt.get('neighbors'));
  var filt = ee.Filter.lt('distance', 1000);
  pt=pt.set('neighsSize1_2003', neighs.filter(filt).size());
  filt = ee.Filter.lt('distance', 2000);
  return(pt);
});

var campcount_3 = pointdata_withbandcounts_3.select(
    ['neighsSize1_2003']
);

print(campcount_3, "camp count 2003")

I have been using insights from here to try to write a function that executes the above function over all eleven years. So far it looks as follows:

// indices for years you want to calculate stats for (2001 to 2014)
var indices2 = ee.List.sequence(2011,2012);

var spatialFilter = ee.Filter.withinDistance({
  distance: 30000,
  leftField: '.geo',
  rightField: '.geo',
  maxError: 10
});

// map over feature collection of tile centroids
var centroids_new = tileCentroids.map(function(feature){

//Within this function, iterated through the 2001-2012 index
var cdata = indices2.iterate(function(i, dictin2){ //Start of second function with ({ open

//Generate the string name of the variable in the camp data (open_2001 = camp open in 2001)
var operatingThisYear = ee.String('open_').cat(ee.Number(i).int());

var joinedUp = ee.Join.saveAll({
  matchesKey: 'neighbors',
  measureKey: 'distance',
  ordering: 'distance'
}).apply({
  primary: tileCentroids,
  secondary: camps.filterMetadata(operatingThisYear, 'equals', 1),
  condition: spatialFilter
});

var withNearestDist= joinedUp.map(function(f) {
  var nearestDist = ee.Feature(ee.List(f.get('neighbors')).get(0))
      .get('distance');
  return f.set('nearestDist', nearestDist);
});

var neighsSize1 = ee.String('camps1km_').cat(ee.Number(i).int());

var getdata = withNearestDist.map(function(pt) {
  var neighs = ee.List(pt.get('neighbors'));
  var filt = ee.Filter.lt('distance', 1000);
  pt=pt.set(neighsSize1, neighs.filter(filt).size()) ;
    return(pt);
});

return ee.Dictionary(dictin2).combine(getdata);

}, ee.Dictionary()   );

 return feature.set(cdata) ;

});

Unfortunately, the function to generate 11*20 variables returns this error: FeatureCollection (Error) Dictionary: Unrecognized arguments to Dictionary constructor. What is the error that's making this function not execute?
Full code here.


texts:
Does anyone know of a way to silently install an ESRI Data License on a machine?

So far, I've tried:

Creating the directory C:\Program Files\ESRI\DataLicense and placing
the *.sdlic file in there as indicated in this legacy forum posting.
Running the DataLicInstall.exe that ships with ArcGIS for Server
against the *.sdlic file, but I get an error "Initialize ArcObjects
Failed".
Using AutoHotKey to manually install the license via
simulated button clicks but unfortunately the Windows UAC button will
not accept simulated button clicks from AHK. So when the UAC is raised the user still
needs to press Yes to allow the ArcGIS Administrator to make changes to the computer.


texts: Use Python to import MXDs into an ArcGIS Pro Project and add the MXD file name to the imported maps

texts:
The method of importing an MXD to an ArcGIS Pro project has already been asked as Importing MXD to ArcGIS Pro programmatically? and it goes like this:

workspaceMXD = r'FOLDER CONTAINING MXDS'
projectsArea = r'C:ArcGIS\Projects'
APRX_name = 'ArcGIS Project.aprx'
output_APRX = os.path.join(projectsArea, APRX_name)

# Variables
arcpy.env.overwriteOutput = True

aprx = arcpy.mp.ArcGISProject("Current")
arcpy.env.workspace = workspaceMXD

# List files in workspaceMXD
listFiles = arcpy.ListFiles("*.mxd")

# Import MXD Files
for files in listFiles:
  # file names in list appended to folder path
  inputMXD = os.path.join(workspaceMXD, files)
  print (inputMXD)
  aprx.importDocument(inputMXD)
  aprx.saveACopy(output_APRX)

This results in a folder of MXD's being imported into an ArcGIS Pro Project - however, these imported maps are named according to their data frame name (e.g. Layers).

How can you access the MXD filename and include it as a variable in the aprx.importDocument class?


texts: How to use FME to parse .txt file to publish content to REST Interface?

texts:
I have a, from my point of view, complex question, which can be hopefully answered by some  of you. Firstly some general information. I'm using

Windows Server 2012
FME Desktop and Server 2013
ArcGIS Server 10.2

My plan: I want to use a Web Harvesting program, which searches on OpenData-Sites for accessible WMS Services, which I want to integrate in my WebMap. This application is already working. It searches for available WMS Services using defined filters and stores the information in a .txt file. So the .txt file stores all the necessary the WMS-Links
(a link could look like this:

http://gateway.hamburg.de/OGCFassade/DE_HH_WMS_INSPIRE_A3_4_Bodennutzung_lokal.aspx?SERVICE=WMS&REQUEST=GetCapabilities&VERSION=1.3.0

or this

http://www.geoportal.rlp.de/mapbender/php/wms.php?REQUEST=GetCapabilities&SERVICE=WMS&VERSION=1.1.1&layer_id=32327

Now I want to use FME Desktop/server, which reads/parsed the lines of the .txt file (the lines are the WMS links) and after that, FME should publish these links to the FME or ArcGIS Server Rest Interface, where they are accessible for using it in a webmap.

So....firstly my question would be: Is this possible?
And secondly: If it is possible, how to realize it?


texts: Set custom spatial reference to base map in arcgis javascript api

texts: How can I set my spatial reference to the base map in arcgis api. I want to change the default spatial reference.


texts:
I'm using QGIS 2.18 (Las Palmas).

I want to change the length of the decimal place, I tried the approach in How to change length to 2 decimal places? however, it still displays 13 decimal places.


texts: What is the size of a "point" unit in Esri software?

texts:
When creating point symbology with Esri software, we can use points, in, cm, or mm as a unit of measurement.

Due to the fact that there are a few definitions of a "pt" or point, which of the following definitions should be used?

Didot point (0.376065 mm)
Pica point (0.3514 mm)
DTP point (0.3528 mm)


texts: Change text\ font size element using Arcpy

texts:
I try to change the font size and the element xy position of text element for several maps. Right now the font size is 21.79 and the xy position is:

elm.elementPositionX = 8.204

elm.elementPositionY = 26.596

I try this code:

import arcpy

mxd = arcpy.mapping.MapDocument(r"G:\desktop\Project\project.mxd")
for elm in arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT"):
    if elm.text == "data1":
        elm.fontSize = 16
        elm.elementPositionX = 9.898
        elm.elementPositionY = 27.649
        mxd.save()
del mxd

The code above works fine when i run for one map but when i loop it on several maps the code does not work:

import arcpy, os, sys
from arcpy import env

env.workspace = r"G:\desktop\Project"
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"G:\desktop\Project\\" + mxdname)
    for elm in arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT"):
        if elm.text == "data1":
            elm.fontSize = 16
            elm.elementPositionX = 9.898
            elm.elementPositionY = 27.649
            print 'elementPosition changed'
        mxd.save()
del mxd


texts: Keep getting NAs when I run sp's over function on LatLon points on shapefile polygons

texts:
I am not sure why i am keep getting NA whenever I run the Over function with Latitude and Longitude point on the polygon from shapefile. Please note that this is first time for me doing the spatial analysis, but I have done my research and replicated things, but didn't succeed. I need some points which are outside of the polygon to be NA, so I can focus on the real data.

I read these sources since these pertain to my cause but I can't work my problem out:
https://stackoverflow.com/questions/19002744/spover-for-point-in-polygon-analysis
Checking if points fall within polygon Shapefile
R: Error in checking point inside polygon

Here is my code chunk

library(sp)
library(rgdal)
library(readr)

gainsville_df <- read_csv("311_Service_Requests__myGNV_.csv")
gnv <- readOGR("~\\Downloads\\GIS_cgbound", layer = "cgbound")

gnv_latlon <- spTransform(gnv, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

gnv_raw <- data.frame(Longitude= gainsville_df$Longitude, Latitude= gainsville_df$Latitude)

coordinates(gnv_raw) <- ~Longitude + Latitude
proj4string(gnv_raw) <- proj4string(gnv)
over(gnv_raw, as(gnv,"SpatialLinesDataFrame"))

#Yeilds:
#  FID_cgboun Id Perimeter Area Acres Hectares Shape_Leng
#1         NA NA        NA   NA    NA       NA         NA

# Desired Output:
# Whereas I should have seen which gainesville Latitudes and Longitude are within the shpaefile
# polygon so I can drop the outliers, that have the NA. According to this, none of my LatLon points
# are inside the polygon.

The datafiles are here:
Shapefile: https://github.com/THsTestingGround/SO_readOGR_quest/tree/master/GIS_cgbound
reading csv file: https://github.com/THsTestingGround/SO_readOGR_quest/blob/master/311_Service_Requests__myGNV_.csv


texts: MapBasic tool not compiling due to undefined variable

texts:
I am attempting to compile a quick MapBasic tool I have written, however it won't compile despite the variable being clearly defined

It's probably something simple I am missing, can somebody point out where I am going wrong?


texts: Averaging Cartesian ECEF coordinates then converting to WGS84 without latitude drift

texts:
Preface: I'm unsure if GIS or math is the correct community for this question.

I have a set of 9 WGS84 latitude/longitude (degree) pairs for which I need to find the weighted average of.

To prevent things like the dateline causing massively warped averages, I'm converting from WGS84 into ECEF Cartesian (X, Y, Z), weighted averaging those coordinates, then converting back to a WGS84 latitude/longitude (degree) pair.

The problem is that this produces an incorrect result - the latitude seems to 'drift' the further from the equator you are.

The conversion algorithms used (from WGS84 to ECEF, then ECEF to WGS84) work as intended and produce accurate results. The problem appears to lie with the weighted averaging of the Cartesian coordinates - I'm not sure why this is.

Question: How do I produce a weighted average of a set of Cartesian coordinates that, when converted back to WGS84, will not cause the latitude to drift?

Notes:

Naively weighted averaging unconverted WGS84 latitude/longitude pairs produces the correct result (unless the set crosses the dateline, for example).
As a test, I converted the naively weighted averaged WGS84 lat/long pairs to and from ECEF, and this produced the correct result, which means it's definitely something to do with weighted averaging the ECEF coordinates and not the conversion algorithms used.

Sample code (in Python):

WGS84 -> ECEF:

lat = np.radians(lat)
lon = np.radians(lon)

d = e * np.sin(lat)
N = a / np.sqrt(1. - d * d)

x = (N + ele) * np.cos(lat) * np.cos(lon)
y = (N + ele) * np.cos(lat) * np.sin(lon)
z = ((1 - esq) * N + ele) * np.sin(lat)

return x, y, z

ECEF -> WGS84:

r = np.sqrt(x * x + y * y)
Esq = a * a - b * b
F = 54 * b * b * z * z
G = r * r + (1 - esq) * z * z - esq * Esq
C = (esq * esq * F * r * r) / np.power(G, 3)
S = cbrt(1 + C + np.sqrt(C * C + 2 * C))
P = F / (3 * np.power((S + 1 / S + 1), 2) * G * G)
Q = np.sqrt(1 + 2 * esq * esq * P)
r_0 = -(P * esq * r) / (1 + Q) + np.sqrt(0.5 * a * a * (1 + 1.0 / Q) - \
    P * (1 - esq) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r)
U = np.sqrt(np.power((r - esq * r_0), 2) + z * z)
V = np.sqrt(np.power((r - esq * r_0), 2) + (1 - esq) * z * z)
Z_0 = b * b * z / (a * V)
h = U * (1 - b * b / (a * V))
lat = np.arctan((z + e1sq * Z_0) / r)
lon = np.arctan2(y, x)
xi = np.sqrt(1 - esq * np.sin(lat) * np.sin(lat))
N = a / xi
ele = (r / np.cos(np.radians(lat))) - N

return np.rad2deg(lat), np.rad2deg(lon), ele

Weighted average:

x = np.average(cart_x, weights)
y = np.average(cart_y, weights)
z = np.average(cart_z, weights)

(where cart_x, cart_y, cart_z are lists of x, y, and z coordinates in same order as their corresponding weights)

Sample data:

Input (WGS84 latitude, WGS84 longitude, Weighting):

53.708794718  -0.876990985207 0.0
56.7088004625 -0.876990985207 0.125
53.9512592773 -0.876990985207 0.125
53.708794718  2.12301475927   0.125
53.708794718  -0.634526425941 0.125
50.7087889735 -0.876990985207 0.125
53.4663301588 -0.876990985207 0.125
53.708794718  -3.87699672969  0.125
53.708794718  -1.11945554447  0.125

...converted to X, Y, Z:

3782945.06909 -57907.7178905 5117626.09215
3508729.39335 -53710.140685  5308169.95693
3761162.0512  -57574.272696  5133553.1923
3780791.31608 140156.069292  5117626.09215
3783156.24975 -41898.5617628 5117626.09215
4046681.01116 -61944.87578   4913093.31984
3804659.61535 -58240.1149506 5101607.59494
3774729.99291 -255812.783436 5117626.09215
3782666.14289 -73915.8369984 5117626.09215

Weighted average of Cartesian coordinates:

3780321.97159 -57867.5646271 5115866.05408

Which produces the incorrect result of:

53.7184377538 -0.876990985207

Expected result of (approximately):

53.708794718  -0.876990985207
Difference: 0.00964303580069 0.0


texts: sa.fill creates NoData raster in standalone script ArcMap 10.0

texts:
I am trying to put together a geoprocessing workflow in Python, and I've hit a stumbling block with sa.fill.  When run as a standalone script, the tool produces a raster composed entirely of NoData values.  If I run the same code in the Python window in ArcMap, a correctly filled raster is produced.  My code snippet is:

outfill = sa.Fill(raster,"")
outfill.save(fill)
print "Sink fill complete."

I've tried copying my entire script into the Python window and adding it to a toolbox, but I get the following error in both cases:

Runtime error <type 'exceptions.EOFError'>: EOF when reading a line

The standalone script does not produce this error.  It actually runs through all the steps without incident, but the output is useless as everything after the Fill step is based on NoData values.

Does anyone have any thoughts on this?


texts:
I have a couple of MapInfo programs that use external EXEs for some functionality (e.g. ogr2ogr for layer conversion and OV2 Tools for satnav data output).

In my programs I would like to have a settings dialog where you set the location of these EXEs and then somehow store these so that they are remembered even after Mapinfo has been closed down.

What's the best way to do this?


texts: How to deal with direction (CW, CCW) of circles, in particular with shapely?

texts:
I am translating an input that describes zones by mixing segments and arcs. Arcs are given by a center, 2 points (beginning/end of arc) and a direction (CW or CCW). In a first prototypes using GDAL, I inspected once how the points were sorted after the call to .buffer() and made the assumption that they will always be in the same direction. Then I moved my code to shapely, and it looks like the direction used is not the same, so I'm back to inspecting and hardcoding the internal direction.

But my guess is that this is not the correct way of dealing with circle direction. Should I check after every calls to .buffer() the direction in which the points are ordered ? If so, is there a well known method to compute the direction ?


texts:
I am trying to split .tif image that i have downloaded from satellite imaginary. The .tif image is 114MB and has lot of data and i want to split the image into further small images to work more on the images.

The code I have return for that:

import os, gdal

in_path = 'C:/Users/naresh/'
input_filename = '20180831_181344_0f49_3B_AnalyticMS_DNN.tif'

out_path = 'C:/Users/naresh/output_folder/'
output_filename = 'tile_'

tile_size_x = 50
tile_size_y = 70

ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
xsize = band.XSize
ysize = band.YSize

for i in range(0, xsize, tile_size_x):
for j in range(0, ysize, tile_size_y):
    com_string = "gdal_translate -of GTIFF -srcwin " + str(i)+ ", " + str(j) + ", " + str(tile_size_x) + ", " + str(tile_size_y) + " " + str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(i) + "_" + str(j) + ".tif"
    os.system(com_string)

While i am trying to execute this code i am getting error in AttributeError: 'NoneType' object has no attribute 'GetRasterBand' in band = ds.GetRasterBand(1) line.


texts: How to calculate the difference between points values on each side of the road at specific locations?

texts:
I am working with a data set of points placed on either side of a given road, containing different attribution (forest, height, etc).

I need to be able to calculate the differences between the attribute on each side of the road at specific locations.

Is there a tool in ArcGIS 10.1 that can perform this operation, or are there any suggestions on algorithms, tools, scripts I could use?


texts:
I installed QGIS 1.7.3  using  OSGeo4W installer.I have python 2.7 and pyqt4-py2.7 installed on my system.My system paths are:

  PATH:C:\OSGeo4W\apps\Python27\Lib\site-packages\PyQt4\uic;C:\Python27\Lib\site-packages\PyQt4;C:\OSGeo4W\apps\qgis\bin;C:\OSGeo4W\bin;

  PYTHONPATH:C:\OSGeo4W\apps\qgis\python

When i try to run my .py files i get following error:

  RuntimeError: the PyQt4.QtCore module is version 1 but the qgis.core
  module requires version -1.

what can be the problem?


texts: i wan to delete the Location of a data set in grass gis 6.4.4 using grass command. For customizing grass for our data set, i need to load different location and mapsets of different regions of country.I seen for mapset removal. To remove location is there any command?


texts:
today I discovered this library, but after trying to convert coordinates from one to other system find out that there is some calculation error.

I googled possible solutions for this problem, but without success. I found this link and followed instructions but it doesn't work to me.

https://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/

I would like to convert coordinate from this zone:
MGI / Balkans zone 6...

to google earth system...

I created this function for testing:

        public static void Check(double  x, double y)
    {
        ICoordinateSystem CoordSystemInput = ProjNet.UnitTests.SRIDReader.GetCSbyID(31276);
        var _testCoordinationSystemInput = CoordinateSystemWktReader.Parse(CoordSystemInput.WKT) as IProjectedCoordinateSystem;

        ICoordinateSystem CoordSystemOutput = ProjNet.UnitTests.SRIDReader.GetCSbyID(32767);
        var _testCoordinationSystemOutput = CoordinateSystemWktReader.Parse(CoordSystemOutput.WKT) as IProjectedCoordinateSystem;

        ICoordinateTransformation trans = new CoordinateTransformationFactory().CreateFromCoordinateSystems(_testCoordinationSystemInput, _testCoordinationSystemOutput);

        double[] fromPoint = { x, y };
        double[] p1 = trans.MathTransform.Transform(fromPoint);

        string Lat = p1[0].ToString();
        string Lon = p1[1].ToString();

        MessageBox.Show(Lat.ToString() + "," + Lon.ToString());
    }

First one system under number 31276 is:

PROJCS["MGI / Balkans zone 6",GEOGCS["MGI",DATUM["Militar_Geographische_Institute",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6312"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4312"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",18],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",6500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","31276"]]

Second one is under number 32767 is:

PROJCS["Popular Visualisation CRS / Mercator",GEOGCS["Popular Visualisation CRS",DATUM["WGS84",SPHEROID["WGS84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7059"]],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9102"]],AXIS["E", EAST], AXIS["N", NORTH], AUTHORITY["EPSG","4055"]],PROJECTION["Mercator"],PARAMETER["semi_minor",6378137],PARAMETER["False_Easting", 0],PARAMETER["False_Northing", 0],PARAMETER["Central_Meridian", 0],PARAMETER["Latitude_of_origin", 0],UNIT["metre", 1, AUTHORITY["EPSG", "9001"]],AXIS["East", EAST], AXIS["North", NORTH],AUTHORITY["EPSG","3785"]]

Here are my testing coordinates and link on correct result:

x = 6503462.895
y = 4922042.835

https://epsg.io/transform#s_srs=3908&t_srs=900913&x=6503462.8950000&y=4922042.8350000

x = 2008028.9823
y = 5534233.6783

But this is what I getting (with mistake):

x = 2008594.2781
y = 5534328.1856


texts: Adding and filling GEOMETRY column in new spatialite

texts:
I have three tables. One has no geometry. The second is POINT type, and shows up correctly on the map. The third is, for the moment, a simple checker table, don't worry about it.

I am able to perform a table join using the following query:

select *
from ignmeta_ages join ignmeta_samples on
(ignmeta_ages.sampleid=ignmeta_samples.sampleid)
where ignmeta_ages.sampleid in (select sampleid from ignmeta_samples) and
ignmeta_ages.sampleid not in (select sampleid from
ignmeta_sample_ages_merge)

This correctly outputs a join of the two. However, the GEOMETRY column is empty of data, and I am completely unable to get it to the point where I can load it as a layer and have the points properly positioned.

What commands should I try for making the data update? X and Y coordinates are found in columns in the second table labelled longitude and latitude.

I have saved a copy of this query to a new table named temp. Running this command makes the GEOMETRY column NULL filled, rather than empty:

update temp set geometry = MakePoint(temp.longitude, temp.latitude, 4326)

I get the feeling this should be the way to go, but it's not working properly.

EDIT-----------------
Some progress (I think). After amending my query to name each column rather than using *, and NOT selecting the old GEOMETRY column, I was able to use

select AddGeometryColumn('temp', 'GEOMETRY', 4326, 'POINT', 'XY')

To add a column of NULL. I then used this:

update temp set GEOMETRY = MakePoint(25, 50, 4326) where ogc_fid=325

And was rewarded with one field that imported properly, with spatial data. My question now turns to 1: How do I replace 25 and 50 with longitude and latitude from the table? 2: How do I have this go down the table and apply it to each row respectively?


texts: Coordinate Conversion - identifying source spatial reference

texts:
I recently acquired a shapefile that when added in ArcMap 10.2, the software says it has an unknown spatial reference system. The X and Y coordinates are both 6 digit coordinates (for example, 712568.799, 257083.471).

I would like to be helped on how to convert these coordinates to Arc 1950 UTM Zone 36s which is a system used in Malawi. If the Y coordinate had 8 in front it could have worked perfectly but this is not the case.


texts: Error in Raster Calculator when copying syntax typed in text editor

texts:
I have a very long operation string using conditionals and simple operations. However, I get an error when copying the syntax to Raster Calculator, apparently the commands are not recognized when copied from text editor.
The problem is I need to calculate about 70 operations with varying values, is there a way to work the syntax outside the command box and make it work?

This is my string and it has been cleaned up, but still get error 00539

Con("slope_dgr",1,0, "VALUE <=20") *
Con("dem_clip",1,0, "VALUE >=1") *
(("acacia_sh_" * 0.6 +
    "e_open_" * 0.9 +
    "e_wood_" * 0.8 +
    "hum_grass_" * 1 +
    "oth_open_" * 0.8 +
    "oth_sh_" * 0.6 +
    "oth_wood_" * 0.8 +
    "trop_e_grass_" * 0.75 +
    "trop_grass_" * 0.5)/9) *
(("yellow_dpx1_" * 1 +
    "red_dpx1_" * 1 +
    "loams1_" * 0.75 +
    "m_earth1_" * 0.75 +
    "sands1_" * 1 +
    "cr_clay1_" * 0.5)/6) *
(("ig_fels_intr" * 1 +
    "ig_fels_volc" * 1 +
    "ig_maf_intr" * 1 +
    "ig_maf_volc" * 1 +
    "ig_volc" * 1 +
    "metamorphic" * 1 +
    "metased_silic" * 1 +
    "regolith" * 1 +
    "sed_carb" * 1 +
    "sed_silic" * 1)/10) *
("wline_mj_" * 0 +
    "wline_mn_" * 0 +
    "warea_np_" * 0 +
    "warea_pn_" * 0) *
Con("wline_mj",1,0,"VALUE =0") *
Con("wline_mn",1,0,"VALUE =0") *
Con("wanp",1,0,"VALUE =0") *
Con("wapn",1,0,"VALUE =0") *
("r1000" * 1 +
    "r800" * 1 +
    "r600" * 1 ) *
("t27" * 1 +
    "t30" * 1 +
    "t33" * 1)


texts: Reclassifying multiple rasters using ModelBuilder or ArcPy using ArcMap10.1 or 10.2

texts:
I have a folder with more than 500 raster (DEM) and I want to reclassify (equal interval method) each of them but number of classes in each case is determined by the number of 100m elevation bands. i.e number of classes =(max elevation- minimum elevation)/100. I have obtained max and minimum elevation of each raster from GetRasterProerties tool.
How do I put this custom number of classes for each raster?


texts: Skipping TopologyException and continuing operation anyway using R

texts:
I'd like to see the result after an hour of working while substracting two shapefiles. This is the code, but don't ask me what in those files. Whatever it is, I'd like to know how to skip those errors and see the result even if it's fake. Those files are huge.

shp1 <- rgdal::readOGR(dsn = paste0("LPDF/SHP/",NAME), layer = layer1, encoding = "utf-8", stringsAsFactors = F)
shp2 <- rgdal::readOGR(dsn = paste0("LPDF/SHP/",NAME), layer = layer2, encoding = "utf-8", stringsAsFactors = F)

Here, I'm trying to use some actions to perform validation as found on some different QAs:

p_ss = sf::st_as_sf(shp1)
p_ss <- sf::st_buffer(p_ss, 0)
p_pol <- sf::st_cast(p_ss, "MULTIPOLYGON")
sb1 <- as(p_pol, "Spatial")

p_ss = sf::st_as_sf(shp2) #shp7
p_ss <- sf::st_buffer(p_ss, 0)
p_pol <- sf::st_cast(p_ss, "MULTIPOLYGON")
sb2 <- as(p_pol, "Spatial")

sb1 <- maptools::unionSpatialPolygons(sb1, IDs = 1:length(shp1), threshold = 1, avoidUnaryUnion = T)
sb2 <- maptools::unionSpatialPolygons(sb2, IDs = 1:length(shp2), threshold = 1, avoidUnaryUnion = T)

So up to this point it worked; without previous buffer the 'maptools' function fails.
But after an hour this fails:

p_dif <- rgeos::gDifference(spgeom1 = sb1, spgeom2 = sb2, byid = F, unaryUnion_if_byid_false = F, drop_lower_td = T, checkValidity = T)

  Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :

  TopologyException: side location conflict at 505082.5 427876.5

Not sure what else could I use in R? byid=T crashes due to memory lack. Maybe I should try to use simple loop and just do it one by one thru the whole data (and day possibly)?

Can you give me some suggestions?

Using suggestions from comments I've tried to make use of automated topology corrections:

Just installed lwgeom using R only and together with other libraries:

p_s1 = sf::st_as_sf(shp1)
sv1 <- lwgeom::st_make_valid(p_s1)
p_pol <- sf::st_cast(sv1, "MULTIPOLYGON")
sb1 <- as(p_pol, "Spatial")

p_s2 = sf::st_as_sf(shp2)
sv2 <- lwgeom::st_make_valid(p_s2)
p_pol <- sf::st_cast(sv2, "MULTIPOLYGON")
sb2 <- as(p_pol, "Spatial")

made these shapefiles valid. But in the substraction:

p_dif <- rgeos::gDifference(spgeom1 = sb1, spgeom2 = sb2, byid = F, unaryUnion_if_byid_false = T, drop_lower_td = T, checkValidity = T)

it only worked if with unary union (TRUE) and byid=FALSE, not sure about dropping. The unary union make it very fast as well. This is enough for me.


texts:
I need some guidance.  I'm using an ESRI Base Map with Web Mercator projection(Wkid - 102110) in an web application using JavaScript. I'm drawing geodesic circles on the Base Map for a given radius and location and getting Map Scale using the ESRI getScale method.
The whole process works fine when im using an address in Florida.

But when i enter an address in Alaska, although the geodesic circle is correct, the scale is completely wrong.

I read a lot of online posts and everyone talks about how to create the geodesic circles which i already had working but not the scale. According to the scale generated by ESRI for the Alaska circle, the distance between the center of the circle and an end point is 19 miles although i intend to draw a circle of only 10 miles radius.
Also interesting to note that the scale on the bottom left of the base map seems to be correct but the one on on the print template is not.
How do i fix this?


texts: Assigning covariate associated to spatial points to bigger set of spatial points?

texts:
I have two data sets with spatial points (in .csv format): data1 with 220 spatial points with latitude and longitude and data2 with 80 spatial points with latitude and longitude. For data2 I have one covariate indicated the genetic origin of each points. Spatial points in both datasets are not exactly the same.

I would like to assign the genetic origin for spatial points in data1. It seems that I need to define around each point in data2 a square (or other) to be able to associate a genetic origin at each points in data1.


texts: HTTPS post to ArcGIS Online feature layer

texts:
In an external system, I want to invoke an HTTPS post that would create records in ArcGIS Online.
Is there a mechanism in ArcGIS Online that can receive HTTPS posts?


texts: Using Feature Dataset Iterator in ModelBuilder?

texts:
I had posted this on the ArcGIS forums, but thought I may get some further input here.

Having issues acquiring the appropriate input and erase features during the iteration process.

Please see below linked ModelBuilder image.  The model is currently not functional as the inner nested model does not currently send each erase feature class to the Erase function. It merely iterates all of them and then sends the last feature iteration class to the outer model. Hence if I run the model on the three files with the wildcard 040112 and three dates with wildcard 100111 (both recursively through folders), I get one correct erased output, but two incorrect ones.

I can't figure out how to send the feature class along for each iteration of the inner nested model. Clearly I am missing something, but I am unsure what at the moment.

I could be barking up the wrong tree as I am relatively new at advanced modelbuilder stuff. For all I know maybe this is more of a python thing, but I am just starting to get into that and it is over my head at the moment.


texts: Error when creating projected coordinate system

texts:
I have written a function to get a Transverse-Mercator IProjectedCoordinateSystem. However when calling the Define function I get an RPC_E_SERVERFAULT (HRESULT: 0x80010105) and no additional information. I am looking at the examples here and here.

Does anyone know what I am doing wrong?

ISpatialReference get_tmerc(double latitude, double longitude)
{
    ISpatialReferenceFactory3 spatialReferenceFactory = new SpatialReferenceEnvironment() as ISpatialReferenceFactory3;

    var wgs84 = spatialReferenceFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);
    var meter_unit = (ILinearUnit)spatialReferenceFactory.CreateUnit((int)esriSRUnitType.esriSRUnit_Meter);
    IProjectionGEN projection = (IProjectionGEN)spatialReferenceFactory.CreateProjection((int)esriSRProjectionType.esriSRProjection_TransverseMercator);

    IParameter[] parameters = new IParameter[5];
    IParameter parameter = spatialReferenceFactory.CreateParameter(100021); //latitude
    parameter.Value = latitude;
    parameters[0] = parameter;

    parameter = spatialReferenceFactory.CreateParameter(100010); //meridian
    parameter.Value = longitude;
    parameters[1] = parameter;

    parameter = spatialReferenceFactory.CreateParameter(100003); //scale factor
    parameter.Value = 1;
    parameters[2] = parameter;

    parameter = spatialReferenceFactory.CreateParameter(100001); //false N
    parameter.Value = 0;
    parameters[3] = parameter;

    parameter = spatialReferenceFactory.CreateParameter(100001); //false E
    parameter.Value = 0;
    parameters[4] = parameter;

    object name = "TransMerc_Name_ASDFSDF";
    object alias = "Custom_Transverse_Mercator";
    object abb_obj = "Trans_Merc";
    object rem_obj = "Transverse Mercator";
    object usage_obj = "";
    object wgs_obj = wgs84;
    object unit_obj = meter_unit;
    object proj_obj = projection;
    object para_obj = parameters;

    IProjectedCoordinateSystem proj_coord = new ProjectedCoordinateSystem() as IProjectedCoordinateSystem;
    IProjectedCoordinateSystemEdit proj_coord_edit = proj_coord as IProjectedCoordinateSystemEdit;
    //***ERROR IS HERE***
    proj_coord_edit.Define(ref name, ref alias, ref abb_obj, ref rem_obj, ref usage_obj,
                           ref wgs_obj, ref unit_obj, ref proj_obj, ref para_obj);
    return proj_coord;
}


texts: Using data from 30-meter DEM map to fill into missing pixels in 8-meter DEM map

texts: I have a DEM map of 8-meter resolution, where some groups of pixels are blank due to the clouds (I guess). Is there a way in ArcGIS or GRASS to fill into those null areas with data from 30-meter DEM map?


texts: ArcGIS Online's new map viewer - sublayers?

texts:
I just opened up ArcGIS online for the first time in a few weeks and I see their new map view is out.  I'm trying to get comfortable with it, but I'm struggling to find the sublayers option for layers that have more than one sublayer.  For example, for these LANDFIRE layers, there are multiple sublayers you can access.  With the classic map viewer, they're easy to navigate:

Under the content section you can clearly see the sublayers and toggle them on and off.
Now in the new UI:

For the life of me I don't see anywhere you can look through the different sublayers.  I must be missing something obvious here?


texts: How to create folders from shapefile attributes using ArcPy/Python?

texts:
I am wondering if anyone has solved this or can help me on my code. I have a Parcel shapefile, and I would like to create folders based on the attribute data (per row) from a field in the shapefile; the field is FID. The folders will eventually be used for relative hyperlinks.

My error comes in the line os.makedirs(os.path.join(PF,row.FID)); I know the rest of the code works when I comment this out - Python returns all the values per row> I want to create new folders from these values.

import os, arcpy

featureclass = "W:\\Shapefile\\Shapefile.shp"
PF = "W:\\12345678\\Project\\Parcel"

rows = arcpy.SearchCursor(featureclass)
row = rows.next()

while row:
    print row.FID
    os.makedirs(os.path.join(PF,row.FID))
    row = rows.next()


texts: Where can I find GIS dataset (.shp, .gdb, etc) on Washington DC lead pipes?

texts:
I keep seeing articles about DC recently publishing a map of pipes in the city. I can only find the interactive version but need it in a .shp, .gdb or some format that can be used in Arc (.txt or .csv with x,y field/a field that I can relate with a census .shp such as census code or zip code, but preferable the former). Here is one article I've seen: http://dcist.com/2016/06/dc_water_map.php

Anyone know where I can find this?


texts: ArcGIS script for mapping lat/long coordinates to boxes/rectangles?

texts:
Working on a project that requires me to take lat/long coordinates and map them as rectanglular boxes on ArcMap 10.1. To be clear, the upper left & lower right corners of each box/rectangle are determined by two sets of lat/long coordinates, that have been compiled as pairs in the following format in an Excel spreadsheet:

Upper Left X/Y | Lower Right X/Y
--------------------------------
(XX.X, -YY.Y)  | (XX.X, -YY.Y)

(so "(XX.X,-YY.Y)" are in single cells...multi-valued entries, yes.)

I need to map them over the contiguous US, so I'll be using NAD 83 as my spatial reference. All the resulting box/rectangle features should be saved into one polygon Shapefile.

My question is, how do I go about automating this whole process with a Python/ArcPy script, or via some other method? I know this should be simple, but I don't know a single thing about Python, and need to replicate this process for hundreds of boxes.


texts: How to get started with mapping the flow of gas pipe network using GIS?

texts: I'm looking for some specifications, whitepapers and any other related material about urban gas pipe networks in particular about geometric network rules that are special to gas flow.


texts: Finding overlay percentage between two polygon shapefiles

texts:
I have two layers that I am currently working with: 1. Parcels/property boundaries and 2. Wetlands

I am trying to determine the parcels that have 0 wetlands and then for the parcels that do have wetlands, what the percentage of wetland coverage is for EACH parcel?

Is there a simple way to do this? Parcels are in orange, wetlands are in the brighter blue.


texts: How to read a csv dataset with multipoint column in sf

texts:
i have a csv with a multipoint WGS84 column:

head(df, n=3)

          CUP                                      the_geom
        <chr>                                         <chr>
1 F11B00000210001                         (16.578441,40.669874)
2 F11B01000360001                         (16.390963,39.752256)
3 F11B04000210001 ((15.229943,40.441221),(15.889952,40.516008))

I'm trying to get the geometries from it, and i think that "the_geom" could be a simple feature column.
So i try to read it with sf package:

df_ st_read(dsn = 'data,csv', geometry_column = the_geom)

but i get this error:

Error in st_sf(x, ..., agr = agr) : no simple features geometry column present

I think i must change something the column strings, but i don't know so much about sf.

There is a solution?


texts: Is a PK required for GeoPackage vector layers?

texts:
According to geopackage.org's "Getting Started With GeoPackage" guide, features are stored in user-defined data tables. It seems like it says that these tables require a primary key:

  Other than the geometry column and a primary key, the schema of a features table is up to the implementer.

(Emphasis mine)

Is that a hard requirement? The underlying SQLite probably doesn't require a primary key, but does the GeoPackage specification mandate each layer to have one?

Motivation:

I'm asking because QGIS 3.4.4 seems to handle this inconsistently. If you intersect two vector layers the resulting layer will often not have any unique attributes (except for maybe the geometry column if the input layers were topologic).

If I try to store the result directly to a GeoPackage layer, features that would lead to UNIQUE constraint violations on the fid column (which corresponds to the fid of the first input layer and is apparently automatically chosen as the new GeoPackage layer's PK) are simply omitted, turning the resulting layer incomplete.
If I store the result in a "temporary layer", the result is complete. fid of the first input layer becomes fid column of the temporary layer and fid of the second input layer becomes fid_2 there.

If I naïvely try to persist it in a GeoPackage layer using "make permanent", that fails due to UNIQUE constraint violations.
If I use "make permanent" but in the dialog remove fid from the "FID" field, I can save to a GeoPackage layer fine. (And loading from it confirms that all the features are still there and that no new unique column for a PK has been added.)

If I start from a layer loaded from a GeoPackage which apparently has a PK (field fid has "Not Null" and "Unique" checked as restrictions in Layer properties > Attribute form) and in remove the fid field in Layer properties > Source Fields, QGIS refuses to save that layer when trying to leave "editing mode" for that layer.

(Not tagging with qgis, qgis-3.4, etc., as the question itself is about the file format and should be independent of what Software is being used.)


texts: Delineating boundary around many small polygons to reduce number of polygons but retain relative coverage using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.3 with an Advanced license.

I ran the slope tool on a DEM and converted the resulting raster into a shapefile with two polygons, one with slopes > 40% grade and the other with slopes <=40% grade. This creates many polygons across the area I am looking at. I am wondering if there is a way to create an outline/boundary that encompasses the small boundaries, but retains the general shape.

In the pictures above the black rectangles are the small polygons I am talking about. In the second image the blue boundary is what I am trying to create. Ultimately I'd like this blue boundary to be a polygon. The goal of this is to reduce the number of polygons and the size of the shapefile. (I have already used the dissolve tool).