texts:
I want to use an URL of Google satellite for WMS and I don't know from where I can have it.

I tried Bing but seemed to me not up to date and this is a problem for me because I'm collecting data from field and I need an updated map to be able to locate my objects via the existing landmarks.

If you have any ideas about how to have it I would be thankful!!


texts: Select by attribute does not work on an indexed field

texts:
I came across an issue where I can no longer run a select by attribute on a field (In my example it is "Students_7") in a shapefile that previously worked. Even when I build the query using the wizard (see screenshot) it does not return any results.

I tried making a new column in my shapefile, with the same datatype, and then used the field calculator to make my new field equal to Students_7. I then ran the same select by attribute query on this new field and it worked.

Does anyone know why the select by attribute function no longer works on my Students_7 field? The only thing I did today, which is when this change happened, is join some files together, I've since gone and made sure to remove all other joins, I can see that "Students_7" is an indexed field, but I'm not sure how or when this happened.


texts: QGIS3.12 vector drop-down menu does not have openstreetmap as an option

texts: QGIS3.12 vector drop-down menu does not have OpenStreetMap as an option. Processing is enabled, but still nothing... I'm new to this platform and have no idea what I should do to fix it.


texts: How to convert from Large RasterLayer to Formal RasterLayer class in R?

texts:
I'm a bit confused about types here. While I'm loading a image:

tif = raster(paste0(path))

I get some "Formal Class RasterLayer" object. But while I change it's matrix of values, for example:

idx <- values(tif) != 200
  values(tif)[idx] <- 255

if becomes "Large RasterLayer" object. And indeed it's large, for file of about 100Mpx it's over 270mb, while that first is about 10mb in memory.

How can I convert this new one into the "Formal" class to save memory?

I tried: tif1 = raster(tif) but object missing coordinates. as.raster(tif) returns raster chr variable.

EDIT

Found a workaround:

writeRaster(tif, "tif.tif")
tif = raster("tif.tif")

It takes much memory for process but at last its freed. Wish to find something faster without saving to disk.


texts: Natural Breaks - Slice vs Reclassify in ArcGIS

texts:
I'm wondering what is the difference between Reclassify tool and Slice tool while using Natural Breaks method.

Below some output from Reclassify tool - 5 class; Natural Breaks (Jenks) method:

Here output from Slice tool - Slice Method: Natural Breaks, 5 zones

Those two raster are very similar but they're not the same. Maybe You can see it better here:

I'll be very glad if anyone could tell me why those two raster are not the same. They both using Natural Breaks method.

Maybe Slice using different algorithm than Jenks?


texts: Using Title Case on Everything Except Last 2 Characters

texts:
I have addresses in my attribute table that I converted to title case from upper case.

  Example: 123 FAKEVILLE ST NW to 123 Fakeville St Nw

How can I make it so that the last 2 characters in all of the addresses in the Address column are in uppercase? In the example above, I need Nw in uppercase.


texts: QGIS: How to speed up raster processing?

texts:
Is there any way to speed up raster processing in QGIS? I have seen discussion about speeding up rendering and querying, but not processing.
I am using SAGA tools to calculate morphometrics and the calculations are taking extremely long. I am maxing out my CPU. Does QGIS parallel process using all available cores by default?


texts: Has MapInfo text objects formatting tags as ArcMap annotations?

texts:
Esri ArcMap allows you to set annotation formatting via formatting tags, eg:

<FNT name="Arial" size="12">Text size = 12</FNT>
<FNT name="Arial" scale="200">Text scale = 200</FNT>

Through formatting tags we can set one part of annotation with one font style and a different font for other parts.

Does MapInfo have such functionality?
If yes, what is it called and what versions support it?


texts:
I am having trouble using the CCN tool.

My data is a DMS 442315 1214598 and I need to be in DD to work with my county data which I know is in the correct format.

I tried adding a new field then calculating it by using the string

left([Lat], 2) &""& right ([Lat], 2)

etc.


texts: Measuring performance difference between Merge and Append in ArcGIS Desktop?

texts: Often I find myself in a situation where I don't care whether or not my tool produces a new feature class, but I do care how long it takes to combine all of my large datasets. Does it take longer to produce a new feature class using the Merge tool instead of the Append tool, or are Merge and Append essentially the same in terms of performance?


texts: Analyzing spatial differences: conceptual and technical questions

texts:
My question is both conceptual and technical, but more conceptual than technical.  It's a little long and multi-pronged, but I figured the best answers would consider everything, so thanks for bearing with me. I list specific questions later.

Right now I'm research the concept of "academic mismatch" in undergraduate education, or when students attend colleges where their academic background is substantially different than the typical student (e.g., a "little fish in a big pond:" a student with low grades/scores compared to his peers, or a "big fish in a little pond:" a student with high grades/scores compared to her peers). While I'm on solid footing with outcomes analysis (do mismatched students do better/worse in school?), I want to explore the mechanisms of attending a school as a mismatch in the first place.  There are a lot of reasons for this, some obvious and well explored (preferential admissions policies), others less obvious and less well explored (where I hope to contribute).

One question I want to explore, and many others (including my dissertation advisor) have found interesting, pertains to where students live and their likelihood of mismatching. The basic hypothesis is that when a student has fewer colleges to choose from within his "local college market," the more likely he is to attend a college that is less selective than he could because it's convenient. Likewise, for non-selective schools, students sometimes attend less academically appropriate colleges because they can get in, even if they're outmatched by their typical peers. There's some evidence out there in local contexts (e.g., Chicago), but very little in a national context, and almost none in rural contexts (which I/we find most fascinating).

This is the basic question I want to answer:

Does the concentration of "good-fit" colleges that are geographically near to a student associate with that students likelihood of attending a college where her academic background matches her peers?

So from this, I need to do the following:

Create a measure of the concentration of good fit colleges. I have methods for deciding if a college is a good fit between student and college, but assessing concentration is more challenging.
Test whether this measure is associated with students "mismatch status"
Create maps demonstrating patterns found

Here are my questions:

My intuition for concentration is to simply count the colleges that are a good fit within a specific distance (or a market defined by a shapefile) for each student and do appropriate analysis.  What better ideas do you have?
Are there geospatial statistical tests that are better than my plan (logistic regression analysis with outcome being match/not match, independent variable being the measure of college concentration?
Software wise, should I bite the bullet and learn ArcGIS?  I have experience as a web developer so my intuition is to run to web-based software, but I'm not sure where to start.  My problem seems like it could be solved without ArcGIS, and if I can learn software I can use in web development, all the better.

As background, I have the following data:

A student's home zip code and a number of variables on their academic background (high school gpa, SAT, etc.)
A data file of all accredited colleges and their academic profiles (which I can assign geo-coordinates)

Everyone is in the USA.

I have time to learn, but am in proposal-writing stage, so insights in that regard are especially helpful.

As background, I have basically zero experience with GIS, but very good statistical analysis, computer programming skills, and data manipulation skills (for a social science researcher). Put another way, I'm comfortable with research design and Stata, R, Python, MySQL, but have never done anything with GIS (besides some web analytics using geolocation with IPs).


texts:
I want to be able to select a rectangular area of the map, but the selection must have a fixed aspect ratio. Let's assume 2:1 like this:

I haven't found any ol.interaction.DragBox option to do this. Currently my only idea is to modify ol.geom.Geometry on boxend event:

interaction.on('boxend', function (evt) {
  var geom = evt.target.getGeometry();
  //TODO: shorten one side of geom to the specified aspect ratio
  var feat = new ol.Feature({
      geometry: geom
  });
  source.addFeature(feat);
});

Are there any other options? (Here's a jsfiddle to play around)


texts: Setting Variables on startup using PyQGIS?

texts:
via startup.py it is possible to execute python code when qgis starts.

Now I wonder, whether it is possible to set variable values, eg. a string variable which holds a path or a dictionary of values, which could be used by several plugins?


texts:
I've got a square origin-destination matrix with 23 elements, I would like to draw each one of these elements as a point around a circle, to achieve a result like the following:

Can you suggest me a fast method to draw points around a circle?

I'm doing the following but doesn't works:

1) set central point

2) create a buffer with segments=99 and distance=50, uncheck "dissolve the result" option

3) convert shape to line, I get 99 lines

4) calculate the line length/23

5) convert line to points, set "insert additional points" to "No", use the result of line length/23 (1316383,260) in "distance" field

What I get is 396 points where the lines touch each other’s:


texts: How to clip Raster using multiple polygons and Python?

texts:
Possible Duplicate:
  clipping a raster by multiple datasets or polygons

I have polygon layer and DEM and I want to write a program using python to cut the raster into multi-raster based on the polygon features.so please anybody can help??

actually I am using ArcGIS 10.1 full license

actually (clipping a raster by multiple datasets or polygons.) does not work for my case. I have wrote this code:.............
import arcpy
from arcpy import env
from arcpy.sa import *

path name variables

polygonLayer = r'....'
output = r'......'
raster = r'.....'

name of ID field

fieldName = "OBJECKTID"
cursor = arcpy.SearchCursor(polygonLayer)
for poly in cursor:
    mask = poly.OBJECTID

create feature layer of current mask polygon

whereClause = '"' + focusField + '"=' + " ' " + mask + "'"
arcpy.MakeFeatureLayer_management(polygonLayer, 'currentMask', whereClause)

but I got an error with whereClause......... so please any help how can I create feature layer fro each row...... so then I can use them to clip the raster..... thanks in advance


texts: How To Properly Implement a Bounding Box for Shapely & Fiona?

texts:
I am trying to spatially analyze two layers within same GDB file using Python / Shapely.

Another user was kind enough to give me a procedure to go by in trying to reduce the number of senseless Shapley "within" calls by doing the following:

  Read the features in your 1st layer into a list like list(gdb)
  For each feat in this list do a list(gdbStructures.items(bbox=fiona.bounds(feat))) to get a list of
  the features in the second dataset that have bounding boxes
  intersecting the bounding box of feat.
  For each element of that result list, make a a.within(b) test as you were doing.

  Shapely's within() can be expensive if your shapes have a lot of
  detail. The approach above will use the spatial indexes in your GDB
  files to make those calculations only when needed.

Below is my code snippet:

with fiona.drivers():
    try:
        if os.path.isdir(sourceDatabaseFile):

            numBLDGwithinPerennialWaters = 0
            listOfPerennialWaters = []
            listOfAllBuildings = []

            with fiona.open(sourceDatabaseFile, 'r', layer="HydrographySrf") as gdb: # Pointer to HydrographySrf layer of source file
                with fiona.open(sourceDatabaseFile, 'r', layer="StructurePnt") as gdbStructures: # Pointer to StructurePnt layer of source file

                    for hydroFeature in gdb:
                        if (('ZI024_HYP' in hydroFeature['properties']) and (hydroFeature['properties']['ZI024_HYP'] > 0)):
                                listOfPerennialWaters.append(hydroFeature)

                    for waterFeature in listOfPerennialWaters:
                        for building in gdbStructures:
                            list(gdbStructures.items(bbox=fiona.bounds(waterFeature)))
                            if shape(waterFeature['geometry']).contains(shape(building['geometry'])):

                                ...

I have tried to do variations of the list command where the bbox is invoked... such as assigning ListofAllBuildings to whatever is returned from the list command or appending ListOfAllBuildings.  No matter what I seem to do I get zero after the looping finishes (presumably storing on the last List command called).

I do know at some iterations of the loop where I call the list command, I do get non-zero lengths.


texts: I am trying calculate % area a raster is inside a polygon. I have the raster file as ASCII grid and polygon enclosing that grid. How can i calculate % area of a raster file inside a polygon in QGIS?


texts: Proximity Analysis for POIs and Highways

texts:
I am trying to check if a POI (Point Of Interest) is within 2 mile distance (crow's flight) of the US Interstate Highways. I have to do this check for 2000+ POIs. Since I am not very familiar with GIS, I did some reading and searching, but could not find a well defined and feasible methodology.

So far I have the following:

A .csv file for the POIs with latitude and longitude data.
A shapefile for the US Interstate Highways, downloaded from the National Weather Service website http://www.nws.noaa.gov/geodata/catalog/transportation/html/interst.htm

What would be the best way to tackle this problem? I am proficient in Python, and I have to use opensource / free software (I don't think I can use ArcGIS).

A few ideas I came up with:

Brute force approach: converting the shapefile to points with a tool and computing the distance for each POI and the points from the shapefile, and taking the minimum distance, then checking if the result is less than 2 miles or not. This would be probably be very resource intensive.
Converting the shapefile into line segments with a tool and similar to the approach above, computing the minimum distance between each POI and the line segments. Again, I think this also would be resource intensive.
Using an opensource GIS software such as QGIS (I am not familiar with it yet), and hope some kind of plug in exists that does the described proximity analysis.


texts: Using Geopackage JS and Openlayers in Angular

texts:
I need to prepare a web app which is going to load a geopackage file and show the layers in OpenLayers.
First of all I was checking the current libs and also some demo information and found a minimal viewer for OpenLayers and GeopackageJS at the next URL:
https://github.com/ahocevar/openlayers-geopackage-demo
I did follow the instructions, cloning the project, installing the devs and executing via:
npm start

It seems to work until trying to open a geopackage file, which shows the next error in the console:
SQL.js
client:52 [WDS] Live Reloading enabled.
geopackage.js:18 Uncaught (in promise) TypeError: Cannot read property 'getTileTables' of undefined
    at geopackage.js:18
    at api.js:78

It's something realted the global geoPackage object, I was trying to get more insight in the JavaScript library looking at the code, but spending the hole day without result.
How can I make it work before trying to implementate something similar with Angular?


texts: OpenLayers 3: Vector Layers not Drawing during animated panning or mouse-drag events

texts:
It appears that features on an ol.layer.Vector are not painted or drawn during panning animations or mouse-drag events (in OpenLayers v.3.1.1) Does anyone have a work-around?

To see the problem, take a look at: http://jsfiddle.net/StephenM/489aveLr/6/

If you zoom in on the rectangle covering Antarctica and use the mouse to pan around, you will see that the rectangle is only painted after the mouse button is release and the pan-animation has completed, whichever happens last.

I found some references to this behaviour in older versions of OpenLayers but it was apparently fixed in version 2, already. Why am I seeing it here? I have tested in Internet Explorer and both Chrome and Firefox and the behaviour is the same.


texts: Viewshed calculation throws "ERROR 010462: Geodesic visibility exception:"

texts:
I'm trying to calculate the viewshed from a given point on the top of a building and I'm using the following line:

VS2 = Viewshed2(DSM_model, 'fl_Opt_Prospect', VS2_agl_raster, 'FREQUENCY', refractivity_coefficient=0.13, surface_offset=0, observer_offset='OFFSETA', outer_radius='RADIUS2', outer_radius_is_3d='3D', horizontal_start_angle='AZIMUTH1', horizontal_end_angle='AZIMUTH2', vertical_upper_angle='VERT1', vertical_lower_angle='VERT2', analysis_method="PERIMETER_SIGHTLINES")

It works fine only for certain values of RADIUS2. If RADIUS2 is 1800 metres for example, it works fine. 1900 metres and above it gives me the following error:

ERROR 010462: Geodesic visibility exception: write failure: C:\Users\OPTIMT~1\AppData\Local\Temp\ArcGISPro11624\_vs_{66AAE8BB-9549-4CF7-9E09-1E8937DE4C84}116243372.tmp


texts: Hyperlink zoom to feature in ArcGIS Online?

texts:
Trying to build a hyperlink that zooms to a specific location on an arcgis online webapp.

so you're on some sort of asset management database. you're looking at the specs of a specific block in a city. you click 'show on map' or 'go to asset' and it takes you to the asset on a map on a different webpage.


texts:
I have a lot of polygons representing Forest land use type. Is it possible to style these polygons in such a way like on this picture?

Example polygon I have now:

So the polygon doesn't have edges, rather curves. Is this effect possible to obtain without redrawing all the features using a circular string and other specialized digitizing tools?


texts: incorrect skewed re-projected Grid lines in TileMill

texts:
I am working on a project in TileMill for the first time.
I am trying to overlay a Graticule of the British national square grid.

I generated a series of LineString for the grid into PostGIS, which are stored as native EPSG:27700 (OSGB36, BNG)

I only have a small area of interest, so it is only 100x100 orthogonal lines

My grid is accurately rendered in QGIS [image:2], which I check by comparing with the "Authorative source" on Bing maps [image:1] (this is a raster of the 1:50000 British ordnance survey map reprojected to EPSG:3857 for web mapping)

When I add this layer to Tilemill image:3], it is all skewed and incorrectly re-projected from its native SRS to 900913 (Tilemill's native CRS)

I can actually see some of this issue in QGIS also, when zoomed out:
[image:4 & image:6]

But when zoomed in they are accurate: [image:5 & image:7]

Unfortunately, The output from TileMill is always incorrect. [image:8]

Please could somebody give me some pointers on how to get my grid true.

(I do not have the reputation to post more than 2 images)


texts: Define spatial reference of layer in ArcGis JavaScript API

texts:
I get an error "Map: Geometry (wkid: 3301) cannot be converted to spatial reference of the map (wkid: 4326)" when I want to add only graphics layer.
I found that spatial reference for a map is defined either by the extent passed to the map constructor or by the first layer added to the map, so I think that problem is that my layer is defined without property about wkid. My question is how to properly define it.
My code is following (it works when I add wms layer first, but I want to add only the "vectorLayer"):

 var map;
    require([
        'esri/map', 'esri/layers/WMSLayer', 'esri/layers/WMSLayerInfo', 'esri/geometry/Extent',
        'dojo/_base/array', 'dojo/dom', 'dojo/dom-construct', 'dojo/parser',
        'dijit/layout/BorderContainer', 'dijit/layout/ContentPane', "dojo/_base/json" , "esri/request", 'dojo/domReady!'
    ], function (Map, WMSLayer, WMSLayerInfo, Extent, array, dom, domConst, parser, dojoJson, esriRequest, request, xhr) {
        parser.parse();
        var spatialReference2 = new esri.SpatialReference({"wkt": 'PROJCS["Estonian Coordinate System of 1997",GEOGCS["EST97",DATUM["D_Estonia_1997",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",59.33333333333334],PARAMETER["standard_parallel_2",58],PARAMETER["latitude_of_origin",57.51755393055556],PARAMETER["central_meridian",24],PARAMETER["false_easting",500000],PARAMETER["false_northing",6375000],UNIT["Meter",1]]'
        });
        var map = new esri.Map("map",
                {
                    spatialReference: spatialReference2,
                    center: new esri.geometry.Point(550000, 6520000, new esri.SpatialReference({ wkid: 3301 })),
                    extent: new esri.geometry.Extent(350000, 6370000, 750000, 6630000, new esri.SpatialReference({"wkid": 3301})),
                    zoom: 4
                });
//        var wmsLayer = new WMSLayer('http://kaart.maaamet.ee/wms/alus?', {
//            resourceInfo: {
//                extent: new Extent(350000, 6370000, 750000, 6630000, {
//                    wkid: 3301
//                }),
//                layerInfos: new WMSLayerInfo({
//                    name: 'MA-ALUS',
//                    title: 'Maaameti WMS MA-ALUS'
//                })
//            },
//            visibleLayers: ['MA-ALUS']
//        });
//        map.addLayers([wmsLayer]);

        var vectorLayer = new esri.layers.GraphicsLayer();
        var layerUrl = 'http://loom-gis.geo.ut.ee:8080/geoserver/ermas/ows?service=WFS&' +
                'version=1.1.0&request=GetFeature&typeName=ermas:testdata_geopnt&' +
                'outputFormat=text/javascript&format_options=callback:loadFeatures&maxFeatures=300&' +
                'srsname=EPSG:3301&EPSG:3301';
        var loadFeatures = function (response) {
            for (var i = 0; i < response.features.length; i++) {
                var geometry = new esri.geometry.Point(response.features[i].geometry.coordinates[0], response.features[i].geometry.coordinates[1], new esri.SpatialReference({ wkid: 3301 }));
                var symbol = new esri.symbol.SimpleMarkerSymbol();
                var graphic = new esri.Graphic(geometry, symbol);
                vectorLayer.add(graphic)
            }
            map.addLayers([vectorLayer]);
        };
        window.loadFeatures = loadFeatures
        esri.request({
            url: layerUrl,
            content: { f: "json" },
            handleAs: "json",
            callbackParamName: "loadFeatures"
        });
    });


texts: GeoJson based on shapefile : large coordinates

texts:
I'm quite new to GeoJson (for use with d3), and in need of some help. When I convert my available shapefile (.shp) with ogr2ogr to a GeoJson-file, the resulting coordinates are very big numbers (in the range of 100000.000).

The shapefile is accompanied by a .prj file which states the following :

  PROJCS["PCS_CarMAGBOG",GEOGCS["GCS_CarMAGBOG",DATUM["CGS_CarMAGBOG",SPHEROID["GRS80 Mod",6380687.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",92334.879],PARAMETER["False_Northing",109320.965],PARAMETER["Central_Meridian",-74.14659167],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",4.680486111],UNIT["Meter",1.0]]

How do I reverse-engineer this in order to get true latitude and longitude ?


texts: Selection Tool selects random multiple points around the one aimed for

texts:
I have a problem selecting points in ArcGis - the selection feature seems to be misbehaving, not allowing me to select point features precisely.

I've imported the XY data of households as a csv file, transformed it into a shapefile and added a boundaries layer.

I am trying to precisely select points within the boundaries, but when I try to do so manually or when I "select by location", for some reason the selection tool is really misbehaving. When I try to select points, it always selects a whole bunch of random points around the ones I actually want, making it impossible to make precise selections of the points.

I uploaded a short video that illustrates this: https://youtu.be/juw0g7cZK3A

I tried to research online but can't seem to find an answer. I suspect it is perhaps related to the fact that the coordinate system is unknown after the CSV XY import, but adding a standard one has not improved things.

Can you perhaps explain what is going on? I am using ArcMap 10.3.1


texts: Imagery Band order for Natural and False color imagery

texts:
I have an image which has had its band ordering rearranged, giving it a weird false color. I wanted to change it back to a natural color, so I used GDAL and through trial and error got the bands re-organised.

But i couldn't find out which bands where meant to be where.

How do i find which bands are first in the array (gdalinfo does not do this, it just says there are 3 bands, red, green, blue)
Are there any good online resources giving different band combinations for different applications?


texts:
How can I unselect a feature in custom hyperlink onclick?

Basically I have select interaction that changes style of a marker and after I click anywhere on the map, then feature is unselected (goes back to default style). This part works just fine.

Now I have hyperlink "button" which closes overlay and should also perform unselection for selected feature. Unfortunately I can't really find a way to do this.

I have tried to attach selected feature to close button and setStyle to null, but no avail. Also tried to set default style on the marker, but then marker won't change when clicked on it again.

Example can be found at
http://jsfiddle.net/c3o6fy7g/5/

Explanation: red marker should change to pink and overlay is supposed to appear. Clicking on "x" closes overlay and should do unselect on the feature, after which marker should change back to green.

"x" button functionality is given at line 18 aka "closer.onclick = function() {"


texts: Add several Spectral Indexes over an image collection in Google Earth Engine

texts:
I want to add several spectral indexes in one go over a Sentinel2 TOA collection.
If I do only with NDVI the script works. If I add 2 types of NDWI I get the following
"Syntax Error: Unespected token (61:4)"
Please see below the script I am using
// First define an area of interest
var lat = 65.64;
var lng = 34.35;
var point = ee.Geometry.Point(lat, lng);
//var aoi = point.buffer(100000); // Create an area (1km buffer around point)
var country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
.filter(ee.Filter.eq('country_co', 'AF'));
var aoi = country;

Map.setCenter(lat, lng, 5); // Center the map on this location, zoom level 10

var start = '2019-03-29'; // initial date of the image collection
var end = '2019-05-05'; //final date of the image collection

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

// reduce NDVI over an image collection
// Function to calculate and add an NDVI band
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

// reduce NDWI (1-2) over an image collection
// Function to calculate and add an NDWI (1-2) bands

var addNDWI1 = function(image) {
  var ndwi1 = image.normalizedDifference(['B8A', 'B11']).rename('NDWI1');
  return image.addBands(ndwi1);
};

var addNDWI2 = function(image) {
  var ndwi2 = image.normalizedDifference(['B8A', 'B12']).rename('NDWI2');
  return image.addBands(ndwi2);
};
// Download the Sentinel-2 imagery collection
var imgs = ee.ImageCollection('COPERNICUS/S2')
    .filterDate(start, end)
    .filterBounds(aoi)
    .map(maskS2clouds)
    .map(addNDVI);
    .map(addNDWI1);
    .map(addNDWI2);
    ```


texts: Exporting 3d map in format manageable by VRML Java viewer?

texts:
I am new to GIS.

I am working on a mapping project for scuba diving sites, and I need to create 3d maps viewable by users on android tablets, iOS, websites, etc...

So far, I achieved the map result I expected using GoldenSoftware "Didger+surfer".

The problem now is to eport the 3d map in a format manageable by a vrml java viewer.

Is there a way to do this using QGIS?


texts: What is causing these streaks in my output tif using arcpy.CopyRaster_management?

texts:
This is the line of code that is creating these ...

arcpy.CopyRaster_management(raster, fname, "DEFAULTS","0","","","ColormapToRGB","8_BIT_UNSIGNED")

and I did have a colormap file with specifics to make it that light blue. I don't know where the other color is coming from or the apparrent "streaking". I think I saw something similar posted about this when I was searching for information about this type of conversion but I can't seem to find that post now. The question remains "Why is this happening and what can I do about it"? It seems like it must be a bug. Note that I'm using 10.1 ArcGIS.

This is the INPUT. Its actually a tif made from the same input GRID by direct "data export" from within ArcMap.

This is the OUTPUT


texts: Error "unsupported operand type(s) for +: 'int' and 'NoneType' for variance in python

texts:
I recently came up with an error in a script I am working with in Python for ArcGIS 10.1. The strange thing about it is that this script works for some files and not for others.
Any ideas as to why this is happening?

The part of the code I get the error in is:

data_elevation=[]

with arcpy.da.SearchCursor(path, "Elevation") as sCursor:
    for row in sCursor:
        data_elevation.append(row[0])

var_elevation = numpy.var(data_elevation)

path in code stands for a FeatureClass (e.g. C:/test_fc.shp) containing field "Elevation" (e.g. 140, 140, 141, 150, ....).

Error:

  Traceback (most recent call last):   File "...", line 90, in
  var_elevation = numpy.var(data_elevation)   File "C:\Python27\ArcGIS10.1\lib\site-packages\numpy\core\fromnumeric.py",
  line 2561, in var
  return _wrapit(a, 'var', axis, dtype, out, ddof)   File "C:\Python27\ArcGIS10.1\lib\site-packages\numpy\core\fromnumeric.py",
  line 37, in _wrapit
  result = getattr(asarray(obj),method)(*args, **kwds) TypeError: unsupported operand type(s) for +: 'int' and 'NoneType'


texts: Route event layer failing to draw all features

texts:
I am encountering some unexpected behaviour in attempting to use linear referencing tools to produce a route event layer for a road network.

I have produced a polyline route shapefile from a previously existing segmented layer without issues, which displays as expected- all lines have been reproduced and consolidated by UID. I'm then running the Make Route Event Layer tool. When setting parameters for this tool, initially there is an error message displayed next to the Event Table Route Identifier Field (Error 00840 - The value is not a Route Measure Event), however this disappears once I have set the From and To Measure Fields. The tool runs without issue, and generates the Route Event Layer. In the attribute table, all looks as it should- the correct number of records, all attributes have been maintained, etc. This is when the irregular behaviours start occurring however:

The layer seems to only display a portion of the features. Panning and zooming around tends to change which features are displayed.
When overlayed on the routes layer, the identify tool will display the attributes from the Route Event Layer even when I click on a section which displays only as the parent route layer
When zoomed to full extent and a rectangle selection is dragged over the route event layer, only the segments that are displaying are highlighted as selected in the attribute table. Reversing the selection and creating a new layer creates a layer with all of the expected features in the attribute table, but no display in the project window.
When I export the event layer, again only a small portion of features will display, even though the attribute table contains all features. Interestingly, it seems to be a different subset of features from the event layer- each time I perform the same operation, it seems that this subset of features chosen to display changes arbitrarily. This also changes the behaviour of the identify tool- clicking on a parent route segment with the exported data layer turned on will only show attributes of the parent route segment. Additionally, the features displayed appear to be fixed, in that panning and zooming will not alter what is drawn. Opening this exported layer in QGIS displays the same set of features, again with a full attribute table.

The only conclusion I can draw from this is that something is occurring in the initial caching when the data is first drawn, both as an event layer and as an exported layer, preventing the full set of features from caching to display correctly. Can anyone shed some light onto what might be happening?

I am using ArcGIS Desktop v10.3

Another user seemed to have a similar issue here: https://community.esri.com/thread/44116

Figure 1:

Two different instances of layers exported from the Route Event layer (Blue and Green)
The original route layer (Yellow).

Expectation would be that both blue and green would overlay all yellow lines. Additionally, blue and green differentiate even between themselves despite having been derived from the same event layer using exactly the same process.

Before export, the Route Event layer looks similar to the blue and green layers, but using the identify tool on any of the yellow lines gives the attributes of the Route Event layer rather than the original route layer that it is being displayed as.

EDIT: Running calculate geometry on the output features shows that a lot of them have 0 length, which doesn't explain why the attributes would be showing when clicking in the area they should be, but does explain why they might not be displaying. I'm not sure why that would be though- there does not seem to be any identifiable pattern of difference in the from and to values between those that do have measurable length and those that don't.


texts:
I have a query where I'm choosing one polygon (parcel) that intersects a point. The problem I'm having is that when I run into a parcel that has multi-geometry, it often chooses the wrong polygon because I'm limiting the query to only return one result. How can I make this query "explode" this parcel and only choose the one polygon that intersects the point?

Here is my query:

SELECT  b.gid,
        st_astext(st_multi((st_dump(st_union(b.geom))).geom)) as area
FROM pointschema.point a
JOIN parcelschema.polygon b on st_within(a.geom, b.geom)
LIMIT 1

here is what the polygon looks like (you see the point in the NE corner):

and here is what the result of the query is (orange):

for reference I'm using:

PostgreSQL 9.5.1
PostGIS 2.2.1
Polygon geom: geom geometry(MultiPolygon,4326)
Point geom: geom geometry(Point,4326)


texts: FME: Join attributes in text file with attributes in shapefile

texts:
I have a shapefile with points and different attributes - one attribute MyID. I also have a plain text file with several attributes - each in a separated column:

MyID MyValue
 1    A
 2    B
 3    C
 :    :

I would like to use FME to append my attributes in the text file using the attribute key/value MyID.

Is this possible, and in that case how?


texts: Map disappearing when zoom in or panning using ArcGIS API for JavaScript 3?

texts:
I am using ArcGIS API for JavaScript 3.19.

I have added a dynamic layer without a base map, without the extent, its working , i can do map activities . But when i added extent in the map while creating , it is showing map, but the moment i started zooming or panning , suddenly map disappears . Totally weird.

 'use strict';
import angular from 'angular';
import uiRouter from 'angular-ui-router';
import routing from './map.routes';
import SwiftMap from '../core/swiftmap';
import config from '../../config/config';

/*@ngInject*/
function mapDirective($rootScope,MapService,MapFactory,$http) {

    return {
        template: require('./mapDirective.html'),
        restrict: 'E',
        transclude: true,
        scope: {
            onmapready: '&'
        },
        link: function(scope, element, attrs) {

             var _extent = this.initialExtent;
           currentMap = new Map("mapDiv",
        {
            logo: false,
            extent: _extent,
            zoom: 13,
            sliderPosition: "bottom-left"
        });

            currentMap.on("load",function(){
                MapService.setData('map',currentMap);
                MapFactory.setCurrentMap(currentMap);
                MapFactory.createMapTask(currentMap);
                $rootScope.$broadcast('map-ready',currentMap);
            });
            currentMap.on("layers-add-result", function (evt) {
                console.log(evt);
            });

            dynamicLayer= swiftMap.addDynamicLayer(MapFactory.currentDynamicMapLayer,true,"dynamiclayer1", layerIds,layerDefs);
            MapService.setData('layer',dynamicLayer);
        }
    };
}

export default  mapDirective;

but when i comment the extent , it is working properly


texts: Dynamic Linear Referencing of events in QGIS from Excel or CSV using virtual layer

texts:
I have a PostGIS table with a geometry type LINESTRING M to represent routes.
Here is the table specification:

CREATE TABLE routes
(
  fid serial PRIMARY KEY,
  r_name text NOT NULL,
  geom geometry(LineStringM,31370)
);

For simplicity in this fictional example, the start measure is 0 and the end measure is 1000 for all routes.

I also have an events.xlsx file with the position of my events:

In QGIS, how can I dynamically locate my events from my excel file without uploading it to PostGIS? (my users have only read access to the database).

I've already imported the Excel file in QGIS (with drag and drop):

I used to use the Event Layer plugin but it does only work for QGIS versions <= 2.8, so I now try to do it with a Virtual Layer (Layer -> Add Layer -> Add Virtual Layer...).

I'm first importing the routes and events layers in the virtual layer dialog, and I then write the following query:

SELECT
    r.fid,
    e."Route",
    r.r_name,
    e."Event",
    e."Measure",
    ST_LocateAlong(r.geometry, e."Measure") geometry /*:pointm:31370*/
FROM
    events e
JOIN
    routes r ON e."Route" = r.r_name;

The layer is created but the geometry is not recognized.
I can verify that the join has correctly been executed because I see that the resulting attributes seem to be correct:

However, there is no geometry. The geometry column is not recognized and is always NULL. The Measures are within the bounds of my LinestringM's start and end value.
I've tried to set the geometry manually in the Virtual Layer dialog, but it does not work.

Does anyone know how to make the ST_LocateAlong function work in a QGIS Virtual Layer?


texts: How to programmatically change an SVG symbol in the Vector Layer?

texts:
I'm trying to change the symbol used in a layer via Python. I've been able to tunnel through a vector layer to the symbol layer, but I don't know how to change the symbol it's using. I'd like to update it using a filename.

Here's what I have so far:

qgis.console.clearConsole()
layer = qgis.utils.iface.activeLayer()

r = layer.rendererV2()
sym = r.symbols()[0]

symLayer = sym.symbolLayer(0)

symLayer seems like it should have what I'm looking for, but I haven't found it yet.

Is this even possible?


texts: export spatial information from table that is joined to a shapefile in QGIS 2.0

texts:
I have a table with data values per district for different years. The table looks like this:
district name / year / value. The same district is thus present in the data multiple times for different years. In QGIS I joined this table to a districts polygon Shapefile. THis creates a table with for every year a row with district name and value. Can I export this table in a way the spatial information of the district is stored in the table?

Thanks in advance!


texts: Get CSV of specific data from NOAA NAM-NMM Dataset using wgrib2

texts: I am trying to get data from the NAM-NMM Dataset and am having a lot of trouble getting the data I want out of the grb2 files. How can I get a CSV file containing information like Probability of precipitation, Probability of freezing precipitation, TMAX, TMIN, etc... plus time of forecast for all latitudes and longitudes.


texts: Storing SRID:4326 data as Geometry or Geography?

texts:
I used OGR2OGR to import SRID 4326 point data into MSSQL Server.  Didn't realize that it defaults to creating a geometry column instead of detecting the data as 4326 and creating a geography column.  I know the code to get around this so no issue there.

The data as is in a geometry column shows up just fine in QGIS.  So do I actually need to put the data in a geography column?

thanks!


texts:
I've got a spatial table in Oracle and a WMS-T in Mapserver, defined in a layer that query that table. The table has several columns and I pass the column to select as a parameter in OpenLayers, so the layer definition in mapserver is the following:

LAYER
    NAME "layer_example"
    CONNECTION "oracle/user@oracleservice/oracle"
    CONNECTIONTYPE ORACLESPATIAL
    DATA "geom FROM  (SELECT E.%type% as value, E.YEAR as YEAR, E.SHAPE_60M as geom FROM EMISSIONS E
                      WHERE E.level=%level% ) USING NONE"

    TYPE POLYGON

    STATUS OFF
    TRANSPARENCY 100

    VALIDATION
        'level' '^[0-9]$'
        'type'       '^[a-zA-Z0-9]+$' #column name
    END

    PROJECTION
            "init=epsg:4326"
    END
    UNITS METERS
    SIZEUNITS PIXELS
    TOLERANCE 0
    TOLERANCEUNITS PIXELS
    CLASSITEM value
    ### LEGEND DEFINITION HERE ###

    METADATA
            "ows_title" "emissions"
        "ows_srs" "EPSG:3035 EPSG:4326 EPSG:900913"
        "ows_timeformat" "YYYY"
        "ows_timeextent" "2000/2011"
        "ows_timeitem" "YEAR"
            "ows_timedefault" "2000"
        "wms_enable_request" "*"
        "wms_metadataurl_format" "text/xml"
        "wms_metadataurl_type" "ISO19115:2003"

    END
END

The problem is the legend! Since the columns in the table has got different classes of values for their legends, how can I define the legend on fly? I tried to define legends in different files and then INCLUDE "legend_%type%.map", but the variable substitution doesn't work for INCLUDE. Actually, it would be great if Mapserver had this feature.

Anyway, some ideas and suggestions? Hope the question is clear!


texts: Is there a public source with an SVG library of symbols for structural geology mapping in QGIS ?


texts: update “more Info” in ArcGIS Online WebMap popup

texts:
I am seeking to update the "more Info" text that appears for URLs in my AGOL WebMap popup. My query is a bit more complex than Change “more Info” in ArcGIS Online popup.

My dataset has a number of "Map" fields that contain links to various types of maps. I would like the name of map to appear in the for each URL instead of "More Info". I received the dataset as an excel with the name of each map displayed as a hyperlink to the map and this is how I would like it to be displayed in the popup (see image below).

I don't think this is possible but I hope I am wrong!

Here is what currently is displayed:

This is what I want to be displayed:


texts: Raster images overlaying each other but shouldn't

texts: When I open a new map document and drag different raster images into the map the images stack on one another rather than form the complete image of a country. When I try zooming in, the multiple images flip back and forth in the same tile. The images are already georeferenced and have been used before without any problem. Any idea how to fix this?


texts: 'gdal_merge' is not recognized as an internal or external command, operable program or batch file in OSGeo4w Shell

texts:
I am trying to batch mosaic a set of rasters using gdal in OSGeo4w shell command line. When I try to use gdal_merge I am getting error as below.

    run o-help for a list of available commands
    C:\>gdal_merge
    'gdal_merge' is not recognized as an internal or external command,
    operable program or batch file.

I tried to see the available libraries in OSGeo. I found gdal_merge is not installed while running
o-help

C:\>o-help
                   -={ OSGeo4W Shell Commands }=-
  avcexport                               avcimport
  cct                                     cs2cs
  curl                                    derb
  dllupdate                               fax2ps
  fax2tiff                                fftw-wisdom
  gdaladdo                                gdalbuildvrt
  gdaldem                                 gdalenhance
  gdalinfo                                gdallocationinfo
  gdalmanage                              gdalserver
  gdalsrsinfo                             gdaltindex
  gdaltransform                           gdalwarp
  gdal_contour                            gdal_grid
  gdal_rasterize                          gdal_translate
  genbrk                                  genccode
  gencfu                                  gencmn
  gencnval                                gendict
  gennorm2                                genrb
  gensprep                                geod
  geotifcp                                gie
  gif2h5                                  gnmanalyse
  gnmmanage                               gpsbabel`

But I am seeing gdal_merge in other system.
Also, I am able to find the gdal_merge.py in the below directory.

  C:\OSGeo4W64\apps\Python37\Scripts

Why is it not showing in my system, and how do I fix this issue?


texts: Iterating in ModelBuilder to Generate Output Files in Source Folder?

texts:
I have a very simple model that has a feature class iterator (recursive) that runs through a directory and sub-folders looking for line files and then projects them to UTM_16N. The only issue with the model is the destination for the output files.

Instead of having all the outputs going to a single destination I want them to be created in the source folder. So basically the model would dump the output in the same folder as the source file.

There are hundreds of folders each containing an individual line file and they all have the same name. The only identifier so to speak is the container folder name (counties in this case).


texts: In CartoDB, best way to have a table with multiple contributors?

texts:
In CartoDB, I want to allow a group of students to collaboratively add rows to a single table. The issue is made more complex because the table is to contain lines, not just points, so Georeferencing is harder. Ideally, students simply draw the points on a shared public map, and the rows show up in one table in my account.

The Export function for a table won't work, as it shares the data as a file, rather than allowing multiple users to add rows to a single file. Duplicate Table and Duplicate Visualization keep the data tied to a single user. Merge tables creates considerable overhead. I dislike having a shared username/password for 20 people, and that might not even work, depending on how Carto manages lockfiles etc.

The answer may simply be "Make 20 maps, upload the KML, and merge it all", but before I go that route, I'd like to know if there is a preferred way to take input from multiple users into a single table?


texts: How to alter transparency by spatial location

texts:
I'm trying to create a map that has a all layers within a subject area having 0% transparency (full opaque), while features outside that area (defined by a polygon shapefile) having a reduced transparency of say 50%. The ideal result is a map that has a focus on a specific area, but also shows the surrounding area. I'm working with ArcMap 10.

Similar to this image:

I tried using the clip feature in the Data Frame properties but its not what I'm looking for.

Thanks!


texts: Calculate velocity x and y components from Lat,Lon and Time

texts:
I'm trying to calculate velocity components in the X and Y direction given two points:

P0 = (Lat0,Lon0,Time0)

and

P1 = (Lat1,Lon1,Time1)

The points are pretty close together so the
components at P0 should be very similar to the velocity components at P1 (not expecting great line weirdness).

I'm familiar with the Haversine method: calculate the Haversine distance, calculate the starting bearing, use trigonometry to get x and y distance and divide by time. That seems a little overkill and computationally expensive.

Is there a simpler way to estimate velocity components without doing all the expensive trig calculations that would be required to calculate distance and bearing?


texts:
I am using geoscript to reproject a layer in my geoserver.

I can use geoserver/catalog to access a layer's projection with no issue.

However, I can't change it. My code is like this:

var catalog = require("geoserver/catalog");
var Process = require("geoscript/process").Process;
var Layer = require("geoscript/layer").Layer;
var Projection = require("geoscript/proj").Projection;

exports.process = new Process({
  inputs: {
    sentence: {
      type: "String"
    }
  },
  outputs: {
    old: {
      type: "String"
    },
    neo: {
      type: "String"
    },
    temp: {
      type: "Boolean"
    }
  },
  run: function() {
    var namespace = catalog.namespaces[0].alias;
    var myPrj = Projection("EPSG:2008");
    var myLayer = catalog.getVectorLayer("test:counties");
    //catalog.getVectorLayer("test:counties").projection = myPrj;
    return {
      old: myLayer.projection,
      neo: myPrj,
      temp: myLayer.temporary
    };
  }
});

The code I commented line never works however wps doesn't tell me any error.

Do you know how can I fix this?


texts: ArcGIS 10.1 crashes during connection with ArcSDE 10.2 in Oracle 11g database

texts: I have been having this issues of ArcGIS (ArcMap and ArcCatalog) crashing on me every time I try to connect to the databases. I have updated both client and server patches for Oracle as suggested by ESRI tech support but nothing has changed. I cannot even click on the Properties without getting ArcGIS to crash. I am using Windows 7 Enterprise SP1 64-bit OS. Has anyone encountered this and have solution for it?


texts: How to deploy services (WMS and WFS) in port 80 and geoserver admin in diferent port(e.g. 8080)

texts:
How can I configure my WMS and WFS to be in different port of my GeoServer admin?

The idea is to have a GeoServer admin internal which provides public services. In other words the idea is to configure my server application allowing internal access to GeoServer admin and public access to WMS/WFS services.

How can I do this? This is something that I have to set on tomcat( server.xml ) or I may change this configuration directly on GeoServer admin?

I think that this could be done by using "Proxy Base URL" on Global Settings for OGC Services. Am I following the right way?


texts:
When I use rotated horline or vertline like @Ian Turton suggested in his answer to question Rendering sld style(buffer) SEVERE null java.lang.NegativeArraySizeException, then I get breaks in the hatch, like on the image below. I can't find a way to get hatch without break. This hatch is rotated to the basis of an object by 45 degrees (so I can't use shape://slash), with different angles it looks even worse.

SLD suggested by @Ian

<PolygonSymbolizer>
<Geometry>
    <ogc:Function name="buffer">
        <ogc:PropertyName>the_geom</ogc:PropertyName>
        <ogc:PropertyName>Buffer</ogc:PropertyName>
    </ogc:Function>
</Geometry>
<Fill>
    <GraphicFill>
        <Graphic>
            <Mark>
                <WellKnownName>shape://horline</WellKnownName>
                <Stroke>
                    <CssParameter name="stroke">
                        <ogc:PropertyName>COLOR</ogc:PropertyName>
                    </CssParameter>
                    <CssParameter name="stroke-width">2.5</CssParameter>
                </Stroke>
            </Mark>
            <Rotation>Nagib</Rotation>
        </Graphic>
    </GraphicFill>
</Fill>


texts:
I am currently using the leaflet geoCSV plugin to bring through CSV data into the map frame.

I am struggling to be able to customise the popup to only show certain attributes of the CSV and different titles. Using the Bankias Example the popup currently uses the script below to pull through all of the attributes using the titles within the CSV.

var popup = '';
    for (var clave in feature.properties) {
        var title = bankias.getPropertyTitle(clave);
        popup += ''+title+'
'+feature.properties[clave]+'

';
    }
    layer.bindPopup(popup);

The issue I have is that I cannot change the titles within the CSV as they are exported from a database so I therefore need to manually change them in the script above.

Using the Bankias example would it be possible to only show the "Localidad" and "Provincia" attributes and change the titles to "locality" and "province" without having to change them in the CSV?


texts:
After a huge help from Blah238.

I found exactly what is the problem.

I wasted few days checking our data, and trying many things.

The problem that I cannot do intersection between two polylines to give a result a new polyline.

It was working on 9.3. So I am not sure if it is not supposed to work , or why it is not working on 10.1

So, I modified the code to this.
Now this code is generating COM exception on my machine.

Blah can you please help if you can

        static void Main(string[] args)
    {
        //ESRI License Initializer generated code.
        m_AOLicenseInitializer.InitializeApplication(new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeStandard, esriLicenseProductCode.esriLicenseProductCodeAdvanced },
        new esriLicenseExtensionCode[] { });

        var geom1 = "LINESTRING (2 1, 7 1)".ToGeometry();
        var geom2 = "LINESTRING (3 1, 6 1)".ToGeometry();
        var targetOperator = (ITopologicalOperator)geom1;
        targetOperator.Simplify();

        IGeometry intersectGeometry = targetOperator.Intersect(geom2, geom2.Dimension);
        Console.WriteLine(overlapped1.ToWellKnownText());

        m_AOLicenseInitializer.ShutdownApplication();

    }

    static IGeometry GetOverlappedGeometry(IGeometry sourceGeometry, IGeometry targetGeometry, esriGeometryDimension dimension)
    {
        IGeometry overlappedGeometry = null;

        IGeometry localSourceGeometry;
        IGeometry localTargetGeometry;

        if (sourceGeometry.GeometryType > targetGeometry.GeometryType)
        {
            localSourceGeometry = targetGeometry;
            localTargetGeometry = sourceGeometry;
        }
        else
        {
            localSourceGeometry = sourceGeometry;
            localTargetGeometry = targetGeometry;
        }

        var targetOperator = (ITopologicalOperator)localTargetGeometry;

        if (targetOperator != null)
        {
            targetOperator.Simplify();

            IGeometry intersectGeometry = targetOperator.Intersect(localSourceGeometry, dimension);

            overlappedGeometry = intersectGeometry;
        }

        return overlappedGeometry;
    }

Please notice that the code is different that Blah code down , that I am using the dimension 1 is the intersection, and not dim0
Because I want a polyline
What the problem with 10.1?

Why I was able to do it on 9.3?


texts: GDAL Open, specify the data type of the input netcdf

texts:
I'm trying to open a NetCDF using gdal declaring the data type of the input to Float32.
The Netcdf datatype is declared in the metadata as "Float", but GDAL reads this information as Float64 and not as Float32. This reading error misinterprets my decimal values causing rounding errors.
float 32 to float64 (wrong):
[0.9160000518, 0.920000052] > 0.92  -->  [False,  True]  WRONG
flaot32 to float32 (correct):
[0.916, 0.92] > 0.92 -->  [False, False]  CORRECT
Is there a way to specify the input data type such as gdal.Open(netcdf, DATATYPE="") ?
I was checking the gdal documentation but I cannot figure out
Regards


texts: Creating new layer from another one and adding new fields with Qgis processing python scripts?

texts:
I'm developing a script in processing for personal use and have some difficulties with the Qgis API. I bought The Pygis Programmers guide . It has been useful but I have had some difficulties in defining processes. I just want to create a script that writes me the same shapefile , add two fields and then calculate the ratio perimeter area.

This is my code:

##testes=group
##Ratio =name
##layer=vector
##output_file=output vector

from PyQt4.QtCore import QVariant
from qgis.core import QGis, QgsFeature, QgsField, QgsFields
from processing.tools.vector import VectorWriter

vectorLayer = processing.getObject(layer)

provider = vectorLayer.dataProvider()

writer = VectorWriter(output, None, provider.fields(),
                  provider.geometryType(), vectorLayer.crs())

features = processing.features(vectorLayer)
for feat in features:
writer.addFeature(feat)

del writer

#add fields
New_layer= processing.getObject(layer)
fields=QgsFields()
fields.append(QgsField('AREA', QVariant.Double,'',20,3))
fields.append(QgsField('PERIMETER', QVariant.Double,'',20,3))
fields.append(QgsField('RATIO', QVariant.Double,'',20,3)) #Perimeter/Area   (RATIO)
writer = processing.VectorWriter(output_file, None, fields.toList(),
                              QGis.WKBMultiPolygon, New_layer.crs())
del writer


texts: How to show portion of country and country itself on map (child and parent)?

texts:
I want to make a map using arcmap. I want my map look like this:

A complete map of a province of a country (full size)
In the right or left corner of the map the whole map of country should be shown and on that it should focus or highlight the province shown on map


texts: Using multiple Openlayers.Control.SelectFeature controls?

texts:
I am working with Openlayers 2. The problem that I am facing right now is that I have multiple editable layers in my application I want to add, delete, modify features as well as Displaying Feature information regarding features. So that I uses one SelectFeature Control for displaying info and one for deleting feature. when I select control from one of them it results conflicting with each other.

How can i remove this problem?


texts: Creating point features with exact coordinates in QGIS?

texts:
How do I create point features with exact (manually entered) coordinates in QGIS?

I get precise GPS coordinates from a survey team which I need to add to a point layer. What I want: -Add point, type in the coordinates and when pressing enter the point is created where it is supposed to be.


texts: Mapbox Studio state borders on low zoom level

texts:
I'm editing a style in Mapbox Studio based on "Streets" and I'm trying to show the state boundaries in zoom level lower than 4.
No matter what style starting at zoom level 1 to 22 I give to the layers (admin-3-4-boundaries and admin-3-4-boundaries-bg), whenever I go from zoom 4 to z3.9, the disappear.

I changed the filter at the "Select data" tabs and the same happens. I also tried creating a new layer with those lines and they show perfectly, but also disappear .

You have any ideas what can it be?


texts: Satellite imagery options for small/medium sized countries?

texts:
I need to create land cover classifications for a couple countries in central/eastern Asia. The largest country is 50,000km^2 and the smallest is 10,000km^2. What are my options for recent imagery? At this stage resolution is not my primary concern, so Landsat would be acceptable, but the scan line errors in Landsat 7 prevent me from collecting useful imagery from 2003 onward.

Likewise, I am willing to pay for imagery, but very high resolution data may require too many images to feasibly process.

My goal is to create a land cover classification that can be compared with a previous classification (~1990) created from a composite of Landsat imagery and aerial photos.

EDIT:

It appears that Landsat 8 data is now available through Glovis - imagery goes back to about April 2013. Finding cloud-free imagery is a bit difficult given how few images there are.


texts: How do I disable Mouse Wheel zoom functionality in ArcGIS Silverlight?

texts:
I have an ArcGIS Silverlight app with a Map object.  It containes an ArcGISTiledMapServiceLayer and a GraphicsLayer. We've written a custom measure tool for adding in a point at an offset on the graphic layer.  I am trying to disable the mouse wheel zoom-in/zoom-out function when this tool is activated to prevent the user from being able to zooom while performing a measurement.

I've played around with various combinations of setting the map.IsHitTestVisible = false, graphicLayer.IsHitTestVisible = false, overriding OnMouseWheel(), all with no sucess.  It either had no effect, or it disabled ALL mouse events, including the mouse down and mouse move event which I need for the measure tool.

I've read thru the documentation on the MouseWheel but I guess I'm not understanding it correctly. I just can't seem to get it to do what I want.  And I can't find any specific examples that do what I want either.

Any suggestions or places I can look to figure this out?

I am using Silverlight 4 and ESRI Silverlight 2.4.


texts:
I am working on a project that is modeling runoff at a county level. The DEM used is lidar derived and has a 5 ft resolution. The general workflow I am flowing is

Sink evaluation
Sink Selection (to give a minimum drainage area)
Fill sinks (using the deranged polygon from the sink evaluation and the 'Use IsSink field' option).

When I run the fill sinks tool I receive the following error:

  System.Runtime.InteropServices.COMException (0x80041098): ERROR 010429: Error in GRID IO: MtVatBldNew: Value range for c:\users\user\appdata\local\temp\arcc79c\ah2014~4\tmp2675 exceeds 100000 .
  ERROR 010429: Error in GRID IO: and number of unique values exceeds 500 ..
  ERROR 010236: Error in building VAT.
  ERROR 010067: Error in executing grid expression.
  at ESRI.ArcGIS.GeoAnalyst.RasterConversionOpClass.ToRasterDataset(IGeoDataset dataset, String rasterFormat, IWorkspace pWorkspace, String name)
  at ESRI.APWR.ApHydro.ArcHydroOp.FillSinks(ApLayers& aplayers, Boolean useIsSink, Boolean fillAll, Double fillLimit, ITrackCancel trackCancel, IGPMessage messages, String& exMessage) in C:\Builds\HydroSolutions\10.1_ArcHydro\Sources\ArcHydroTools\src\ApHydro\ArcHydroOp.vb:line 57592
  Failed to execute (FillSinks).

I thought this was telling me that the DEM was too large in size, so I tried running the tool again with several subsets by using a polygon to determine the processing extent. This results in the same error and I am not sure why. I have tried the usual fixes (deleting temp files, checking output locations, using the toolbox vs. the toolbar) and the error still occurs.

In order to get the DEM to a smaller size, it was converted from float to integer so I thought the VAT error pertained to that. I reran the tools using the float based DEM and received the same error.

I have also tried using small subsets of the raster to see if that is the problem, and I have the same error.

Can anyone offer insight on this or help me figure out what I am doing wrong?


texts: Performing meaningful kriging in R (averaging nearby points)

texts:
I have a set of data on the quantitative composition of the particle size distribution of the sedimentary layer. I have to build a map of the distribution of certain fractions across the entire study area.

I want to perform a kriging interpolation over the studied area. However, I have a question about the initial processing of the data set. I attach the map below with the distribution of the sampling points. I wonder if it is necessary to average some points that are very close to each other (especially in places with very dense sampling density)?

Just in case I attach the results of the kriging without any averaging.


texts: How to approximate "Autodesk Geometry Format" string to GeoJson

texts:
I'm having a problem reading following geometry format:

CURVEPOLYGON ((x1 y1 (LINESTRINGSEGMENT (x2 y2, x3 y3, x4 y4, x5 y5), CIRCULARARCSEGMENT (x6 y6, x1 y1))))

...from what I searched this should be AGF Text format used by Autodesk.

My goal is to convert such geometry into GeoJSON (preferably in python3). Since GeoJSON doesn't support curves then approximation is the way to go. Approximation WKT to GeoJSON can be solved.

Does anyone know how to convert such structure into WKT or directly to GeoJSON?


texts: Adding shapefile layers in QGIS causes previous layers to disappear

texts:
This is my first time having such trouble, but I'm only a casual user.

I have a base raster layer of the school I work for and I've been adding vector shapefiles on top of that in separate installments. I made a map of trails and the shapefiles from that show up on the view and layer list. I then add a vector layer of flower beds and as soon as I do that the trails disappear from the main view. The CRS is the same, the Trails are still listed in the Layers and I've brought them to the top of the layer list. I can select them from the attribute table but they still don't show up. Removing the flower beds layer doesn't fix the issue, I need to start a new project to get the Trail layer to show again.

The flower bed layer doesn't cause the topography layer to disappear which is another line shapefile I've reused from another project. I've also added another line layer and point layer on top of the flower beds and they don't disappear.

I've tried this in QGIS Valmiera and Brighton on a Windows XP and Ubuntu machine respectively and the issue is the same.

Any ideas what's going on?


texts: Finding Topology Errors toolbox or extension?

texts:
I do not remember if I add this Topology Errors extension or toolbox to my ArcMap.

I  would like to know more about this toolbox.

This is not the same as the Topology toolbox that comes with the ArcMap


texts: I have I field "MEMORIAL" in one of my vector layers with the source of a html file. The file is a table that I would like to appear in my atlas through HTML frame. How can I do that?


texts: Why does Merge then Dissolve return ERROR 000732: Input Features: Dataset ... does not exist or is not supported?

texts:
I'm not sure why I'm running into this problem. I thought it might have something to do with the merged file being locked or open, but I don't see a "LOCK" file being created in Windows Explorer.

I'm simply trying to Merge a bunch of files and then Dissolve the merged file. Here's my script:

workspaces = arcpy.ListWorkspaces("*", "Folder")
for workspace in workspaces:

    arcpy.env.workspace = workspace
    shps = arcpy.ListFeatureClasses("*", "ALL")

    ## Process: Merge Layers
    descWorkspace = arcpy.Describe(workspace)
    workspaceName = descWorkspace.baseName

    mergeOutput = workspaceName+"_merge"
    arcpy.Merge_management(shps, mergeOutput)

    ## Process: Dissolving Features              <--HERE IS WHERE I'M HAVING PROBLEMS
    disInput = mergeOutput
    disOutput = mergeOutput+"_dis"
    arcpy.Dissolve_management(disInput, disOutput)

This is the error that I'm getting after the merge runs successfully: ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Features: Dataset Bus04_merge does not exist or is not supported

I'm obviously missing a step here... (or using incorrect parameters, but I can't figure out in what way)


texts: Setting extent on web map to pass to geoprocessing service?

texts:
I find example to set map extent that works ok in ArcMap via
..

mxd = arcpy.mapping.MapDocument("CURRENT")
..
df = arcpy.mapping.ListDataFrames(mxd)[0]
newExtent = df.extent
newExtent.XMin, newExtent.YMin = newExtent.XMin+100, newExtent.YMin
newExtent.XMax, newExtent.YMax = newExtent.XMax+1000.0, newExtent.YMax
df.extent = newExtent
arcpy.RefreshActiveView()

I would like to publish this type of  tool to make it generic geoprocessing tool for WebApp Builder- it should take the given web map service and set extent of map in the web session.

It is not clear to me - whether this is possible at all?
How to get web map object with python set its extent ?
If I change the extend of mxd document like to code above, I have impression that all users of the service will get the changed extend, that this is not limited to that particualr web session. In Java script API there is an object Map to manipulate web map, does something like this exist in Python?

I would like to pass map service name as parameter

can pass this as a tool parameter). Is this possible at all?


texts: Proportionally distributing attributes from one layer to another

texts:
I have some experience with ArcMap but only on a relatively superficial level. Right now, I'm trying to perform an action and I'm not sure to get it done. I've tried Googling and documentation but I'm just not exactly sure what function I'm looking for.

I have two layers. One is a grid of equally sized squares, the other is a layer with administrative units. These administrative units all have a value for a particular feature (say no. toilets per household).

Now, I would like to record this feature for the grid layer. When a cell is completely within one of the administrative units, this is straightforward. However, I am not sure how to perform this action when multiple administrative units are contained in a grid cell.

I would like to perform an action that proportionally 'transfers' these values from the administrative units to the grid cells. E.g., say 60% of a grid cell is in administrative unit A, and 40% is in B, with A having 0.5 toilets per household and B having 0.6. I would then like this value for the grid cell to be 60%*0.5 + 40%*0.6 (0.54).

How would I be able to do this?

I have tried spatial join but this does not seem to be able to deal with the proportional nature of my problem.


texts: How to obtain Census Tract Centroids when Bureau's web site unavailable?

texts: Does anyone have the Census Bureau file of census tract centroid, for either 2000 or 2010 tracts? With the Bureau's web site shut down, the only available file is the EPA's, and it omits about 500 tracts.


texts: Combining layers and changing fieldnames in Arcgis Pro ModelBuilder

texts:
I have layers containing polygons. The difference between the layers are the periods of the year about which the layers have information. So layer 1 is called Lawn_p1 and layer 2 is called Lawn_p2 (p stands for period) and so on. The layers contain the same fields. With model builder I want to combine all these layers into one shapefile and also add for every field p1 up to p10. This way for every polygon I can see all the information for every period in the same row.

As a first step I want to change the field names of each layer by adding the period to it: so for layer Lawn_p1 every field should get _p1 at the end. This should then be done for all the other layers as well.
Then I should merge/join the layers.

Does anyone know how to do the first step with model builder for all the layers?

I guess I have to use iteration/for loops.


texts: Selecting points within buffer using PostGIS?

texts: I have approx 73,000 locations. For each location, I need to know which of the other locations is within a 1km radius from it. I would like this result to be in a table (so perhaps using a PostGIS spatial database). This is so that I can create edges between these all locations that are within 1km of each other as part of a network analysis study.


texts: Find all regions (lat, long, radius) that contain a point

texts:
A Region has a latitude, longitude (GEOMETRY('POINT')) and a radius.

Using postgis, how can I query all Regions where the distance from pointParam and region.point is less than or equal to region.radius?

SELECT *
FROM region
ORDER BY region.origin <-> my_location
WHERE (region.origin <-> my_location) <= region.radius // <--- This is where I don't know what to put
LIMIT 5;


texts: Attribute rules for location description

texts:
Our organization is working on migrating to ArcGIS Pro and the new Utility Network. I am working on replicating Attribute Assistant functions with the new Attribute Rules in a test environment.
We have a location description field in all of our utility data that populates the nearest address when a feature is created. I am trying to replicate this function with Attribute rules. Attribute assistant based it on an address centerline file. I am trying to base it off situs points for my first try at the automation.
I am able to grab the nearest address but only one field. For some reason it wont concatenate. If I try to get is to return more than one field I get an error. I have tried concatenate in both the variable and results and about 50 other ridiculous ways to get it to work.
Here is the code that gives me the nearest "street" value. If I try to call and concatenate any other fields I get an error.
So I need the several fields from the intersected layer concatenated together to have a coherent address for my LOCDESC field
var situs_address = "street_name";

var feature_field = $feature.LOCDESC;

var address = featuresetbyname($datastore, 'situs');
//set some variables, tried to concatenate in var = situs_address, arcade does not like that//

var closest_address = first(intersects( address, Buffer($feature, 500, 'feet')));
// this is the meat and potatoes of the operation//

if (closest_address == null)
{
return $feature.locdesc;
}
if (isempty(closest_address.street_name))
{
return $feature[feature_field];
}
//takes care of the null values//

return closest_address[situs_address];
// tried all sorts of combinations here, tried adding variables for each feature i need and concatenating//


texts: Grid options are missing in SAGA toolbox. If I go to options in processing the following message is being displayed. How should I change the path?


texts:
I try to work out a way to assemble the layers I have got in my QGIS project so that no layer appears to be stretched.

original coordinate systems:

vector layer: WGS 84 UTM Zone 34N

raster layer: WGS 1984

In the first case the raster layer appears stretched.

coordinate systems in picture one:

vektor layer: WGS 84 UTM Zone 34N

raster layer: WGS 1984

coordinate reference system: WGS 84 UTM Zone 34N

I can change the coordinate reference system to WGS 1984. What happens is that the raster layer looks fine but the vector layer is stretched.

This does make sense as I have different coordinate systems and QGIS tries to somehow add them together. The problem is that once I try to project any of the layers coordinate systems to the same coordinate systems as the other they appear at very different places.

Is there any possibility to add all layers together without one being stretched?

It is also quite important for me to know if there are solutions to this in general or if there is nothing you can do about it sometimes.


texts:
Using this chart, enhanced my presentations and presentation skills:
 found on Extreme Presentation Method.

How would you modify this chart to be directly related to map design?


texts: How can I stack several 3D elevation layers?

texts:
So far I have been using qgis2threejs to display different layers at once like one terrain elevation layer and one population data layer. The result through is that the 3D population data is not sitting on "the mountains" but on flat   floor with z=0, sometimes hidden below the terrain surface. The layers are displayed in parallel but not on top of each other which should be case.
I could add both layers and display them as one but then I can't manipulation the colors individually with is essential for my work.
I am using qgis2threejs because it is intended for online users to explore. But I am also willing to learn other tools. I am just unsure which can do what I am looking for.
Maybe Mapbox? Deck GL? gVSIG? Do need the commercial ArcScene?
It would really be good if it's possible to explore the resulting 3D visualization on a website.
Below you see a screenshot of my issue. The peaks of population data appear under the mountainous surface, instead of on top of it. I projected the colors on the elevation layer so but thats missing the 3D effect of course.


texts: How to create a shapefile in QGIS using the spreadsheet?

texts: I have a spreadsheet which has columns (X, Y, Att1, Att2 ...). I want to create a point shapefile, all columns in the spreadsheet to be the column in the shapefile. How can I do this in QGIS?


texts: Updating a Field in the Current Selection

texts: I have some address points selected for which that I want to update the address field based upon the row in the selection. How do I figure out the name of the selection taable to pass to arcpy.UpdateCursor()?


texts: Estimate the spatial relationship between two sets of points

texts:
I have two sets of point data, structured like below.

Table A
long | lat | class
 ... | ... | A
 ... | ... | A
 ... | ... | A

Table B
long | lat | class
 ... | ... | B
 ... | ... | B
 ... | ... | B

Both tables aim to capture locations of the same certain things, but Table B is much more granular than Table A. The result is that Table B is quite large (higher n) and Table A is more manageable.

I want to see if Table A is representative of Table B (in other words, is Table A "good enough" to represent table B?). How might I approach this? My thought was to see if they two classes are spatially correlated for a certain geography, but I don't know how I would approach this. I'm using R.


texts: A GeoServer layer is based on PostGis view, but featureID is random (PK is ignored)

texts:
I have a layer in GeoServer.
It is based on PostGis view.
It seems that GeoServer generates random featureIDs.
The view has the same columns as the underlying table, including the PK column.
But GeoServer ignores the PK column for layers based on the view.

When I created a layers based on the table itself, Geoserver used PK as featureID.

How to make GeoServer to do the same for a view?

Thank you in advance.

UPDATE

I've accepted iant's answer. SQL I used to create the metadata table is below. After I ran it on the PostGis DB, I've restarted GeoServer and recreated the layer. It worked.

-- http://docs.geoserver.org/stable/en/user/data/database/primarykey.html#metadata-table-description

CREATE TABLE gt_pk_metadata_table (
  table_schema VARCHAR(32) NOT NULL,
  table_name VARCHAR(32) NOT NULL,
  pk_column VARCHAR(32) NOT NULL,
  pk_column_idx INTEGER,
  pk_policy VARCHAR(32),
  pk_sequence VARCHAR(64),
  unique (table_schema, table_name, pk_column),
  check (pk_policy in ('sequence', 'assigned', 'autoincrement'))
);

INSERT INTO gt_pk_metadata_table(
            table_schema, table_name, pk_column, pk_column_idx, pk_policy,
            pk_sequence)
    VALUES (
    'public',
    'llsoa',
    'identifier',
    1,
    'assigned',
    null);


texts: Transparency not working right in GeoServer layer group

texts:
I have two image mosaic layers of geotiffs. One contains aerial photos, the other contains rasterized vector data, meant as an overlay.

Creating a layer group of these two produces an unpleasant result - jagged edges with black patches:

Expected result (acquired by gluing the overlay to the photo in GIMP):

The overlay is mostly transparent, and the vectors have been rendered using anti-aliasing on alpha channel. When I open just the overlay layer in the built-in GeoWebCache demo, the result is as I expect:

Opening the source TIFF in e.g. GIMP confirms that the anti-alias shading is indeed done against alpha channel, and that there is no background color in the TIFF file.

Close-up of rendered layer group - transparent pixels are not rendered correctly:

What could be causing those black/dark pixels and how might I fix it?

What I have tried:

enabling JAI-ext
setting layers' MergeBehavior to STACK vs. FLAT
playing with InputTransparentColor & OutputTransparentColor

Style is default raster.

EDIT: added actual/expected end results

Output of gdalinfo for one overlay geotiff:

Driver: GTiff/GeoTIFF
Files: overlay.tif
Size is 12000, 12000
... clipped CRS info
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
... clipped coordinates
Band 1 Block=12000x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=12000x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=12000x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=12000x1 Type=Byte, ColorInterp=Alpha

gdalinfo for an aerial photo:

Driver: GTiff/GeoTIFF
Files: aero.tif
Size is 12000, 12000
Coordinate System is:
... clipped CRS info
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  SOURCE_COLOR_SPACE=YCbCr
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
... clipped coordinates
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 6000x6000, 3000x3000, 1500x1500, 750x750, 375x375
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 6000x6000, 3000x3000, 1500x1500, 750x750, 375x375
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 6000x6000, 3000x3000, 1500x1500, 750x750, 375x375

Version information:

GeoServer Version
2.12.0
Git Revision
bc18e09993a30c93f37b7e33e7c54c32b17e014b
Build Date
12-Oct-2017 07:54
GeoTools Version
18.0 (rev 91315ffe07add4d70567a36802e5813e9cc75fa0)
GeoWebCache Version
1.12.0 (rev 1.12.x/6fa967943b81ffff4b46b2ac8d5529f82ed15bed)


texts: ArcGIS Online Map using URL does not work with PhoneGap

texts:
I wanted to show my public accessible ArcGIS Online map with URL on a simple Hello World PhoneGap project with Eclipse. Eclipse threw out the following error in LogCat.

Web Console:RequestError: Unable to load file: //www.arcgis.com/sharing/rest/content/items/{MapId}?f=json status: 0 at http://serverapi.arcgisonline.com/js...3.4compact/:15

Just wondering anyone have had any successful experience with using ArcGIS Online web map with PhoneGap?


texts: Non-spatial SQL Server table with X,Y fields directly in QGIS?

texts:
In QGIS 3.6 I linked a SQL Server table with a field 'X' and a field 'Y', together representing positions in a known co-ordinate system.

How can I show these positions directly in QGIS, without having to use the tool 'Create points layer from table' in the Processing Toolbox (under 'Vector creation') and then using this points layer to view the data in the map?


texts: Splitting LAS file into smaller LAS files by shapefile geometry boundaries - using PDAL or Global Mapper.

texts:
I've looked up pdal documentation on how to clip LAS point cloud files:-
https://pdal.io/workshop/exercises/analysis/clipping/clipping.html

Referencing this example, and having a quick look around, it seems like there only is an option for one output file. Please correct me if I'm wrong.

I am after the input LAS files being split into four separate output LAS files (following the geometries of the indextile shapefile).

I was able to achieve this with the LASTools command lasclip.exe, however I've reached the limit of the 'unlicenced' version of lasclip:-

lasclip -i C:/a/test/BY17_1401merged_test.las -split -poly C:/a/cropped_topo50.shp -o newfile.las

https://www.cs.unc.edu/~isenburg/lastools/download/lasclip_README.txt

In the stdout for the above command;
WARNING: unlicensed. over 3 million points. output slightly distorted.
         tiny xyz noise. points permuted. intensity, gps_time, user_data
         and point_source_ID zeroed.

In the lasinfo report for one of the output LAS files.
WARNING: real max x larger than header max x by 0.020000
WARNING: real min x smaller than header min x by 0.020000
WARNING: real max y larger than header max y by 0.020000
WARNING: real min y smaller than header min y by 0.020000
WARNING: real min z smaller than header min z by 0.010000

Is there any other way of splitting a LAS file into multiple output files using shapefile geometries?
I have a Global Mapper license but cannot find a way to do it with this software. Unfortunately I was unable to get SAGA to work on my laptop and Whitebox is almost unusable (very slow) on my ubuntu computer.


texts: How to select/split polylines in specific distances in bulk?

texts:
I am trying to map road rehabilitation projects. I have polylines (roads) and I have a table which indicates in each year which part of the road was rehabilitated (for example, for a road HW107 rehabilitated from km 75 to km 87 in 2009, km 22 to km 37 in 2010 and etc. I have these numbers in separate columns).

I know how to select/split each individual polyline manually or by indicating specific distance. However, since I have more than thousand road works I was wondering if there is a tool to do it in bulk?

I am using ArcGIS Desktop 10.3.


texts: I uninstalled the ArcGIS 10 from my computer, the software is uninstalled , I failed to uninstall the License manager. It was giving me an error "you don't have sufficient access to uninstall ArcGIS license manager 10" while I logged in as Admin. After reading here, I downloaded REVO software, and I managed to remove it. Surprisingly when I tried to install License manager version 10.1, I got the error  "Another version of this product is already installed, Installation of this version cannot continue". Therefore I cannot install the ArcGIS 10.1.


texts: Using WFS to serve arcs / curves from ArcGIS Server to QGIS?

texts:
I loaded a wfs served by ArcGIS 10.0sp3 into QGIS.
I got complaints that the arc were served as stroked arcs.

Question: is WFS able to serve real arcs instead of stroked ones, or is the client (in this case Qgis or Microstation) responsible for handling the request.? And if so how do I do that?


texts: How to display geotagged photos from my domain in leaflet

texts:
I already asked a question similar to this then I used a plugin of displaying geotagged photos in leaflet from (Google+ photos|Instagram / CartoDB|Norvegiana ) photos in leaflet
This plugin gets photos over the internet connection (Google+ photos|Instagram / CartoDB|Norvegiana ) I want the same but from photos hosted in my domain .


texts: Duplicate/triple a geometry based on a numeric field in the attribute table

texts:
i'd like to duplicate/triple a geometry based on a numeric field in the attribute table.
For example, if there is number 2 corresponding to a record, needs to duplicate this gemometry; if there is number 3, needs to triple this geometry and so on..

Any suggestions, tools or plug-in to do this?


texts: Print large land area using QGIS 2.6.1 Brighton

texts:
I have a rectangular area of approximately 10km x 20km and I wish to print this to A0.
(I will select a scaled portion of this area - most of it - that matches the A0 page size)

I have a few layers in QGIS with OpenCycleMap in the background (lowest layer).

I wish to print a good quality print.

What is the best way to get my features out of QGIS and into a file format that will preserve detail for A0 print while showing crisp definition ?


texts: Merge small underlying polygons into overlaying shapefile with larger polygons

texts:
I have a dataset with the average rent prices saved in a vector layer (.shp). It is "street" based, so it has quite a fine resolution. I have a different dataset (containing info about crime rates) in more coarse (district) based resolution.

I want to achieve that the rent prices are aggregated in the more coarse resolution - say 3 polygons of the rent prices are together one district. I want the mean average rent price of that district.

Anyone got an idea how to solve this?

I'm using QGIS, but R is also fine


texts: Stand alone Calculate Field script using Code Block for IMG Raster in Python

texts:
I am somewhat familiar with using the python language in the Calculate Field tool in ArcGIS Pro, but was wondering how this is done in a stand alone Python script with a raster .img (Combine20012004FltPths.img) file.
I am trying to reclassify a new field I have created using integer values from another field. However, when I try to run the code in IDLE, I get a "SyntaxError: invalid syntax".
Here is my python script for the Code Block and the arcpy management Calculate Field tool:
# Code Block for Calculate Field
codeBlock = def reclass(FlightPat1):
                if (FlightPat1 == 0):
                    return "Unclassified"
                elif (FlightPat1 == 11):
                    return "Open Water"
                elif (FlightPat1 == 12):
                    return "Perennial Snow/Ice"
                elif (FlightPat1 == 21):
                    return "Developed, Open Space"
                elif (FlightPat1 == 22):
                    return "Developed, Low Intensity"
                elif (FlightPat1 == 23):
                    return "Developed, Medium Intensity"
                elif (FlightPat1 == 24):
                    return "Developed, High Intensity"
                elif (FlightPat1 == 31):
                    return "Barren Land"
                elif (FlightPat1 == 41):
                    return "Deciduous Forest"
                elif (FlightPat1 == 42):
                    return "Evergreen Forest"
                elif (FlightPat1 == 43):
                    return "Mixed Forest"
                elif (FlightPat1 == 52):
                    return "Shrub/Scrub"
                elif (FlightPat1 == 71):
                    return "Herbaceuous"
                elif (FlightPat1 == 81):
                    return "Hay/Pasture"
                elif (FlightPat1 == 82):
                    return "Cultivated Crops"
                elif (FlightPat1 == 90):
                    return "Woody Wetlands"
                elif (FlightPat1 == 95):
                    return "Emergent Herbaceuous Wetlands"

#Reclassify field based off another field
arcpy.CalculateField_management("Combine20012004FltPths.img", "NLCD01Clss",
                                "Reclass(!FlightPat1!)", expression_type = "PYTHON3",
                                code_block = codeBlock, field_type = "TEXT")

Would I need to convert the attribute table of my .img (Combine20012004FltPths.img) file into a dBASE or geodatabase table or some other table format, then run the Code Block and Calculate Field tool for it to work?


texts: Snap points to nearest polylines within buffer (fast)

texts:
This is a working version of snapping any point to nearest polyline within a buffer distance. But, this takes significantly longer as each time it copies the point to a different shape and snaps it. Is there a way to create a faster version of this tool?

    for id in ids: #ids contain the list of links within the buffer
        where_clause = """ "_ID_" = %d""" % id
        arcpy.SelectLayerByAttribute_management(base_f, "NEW_SELECTION", where_clause) #selects one/any link
        arcpy.CopyFeatures_management(other_pt_f, m1) #copies it to temp(m1)
        arcpy.Snap_edit(m1, [[base_f, "EDGE", "1 Miles"]]) #snaps the node to that link
        with arcpy.da.SearchCursor(m1, ["SHAPE@XY"]) as curs:
            for xy in curs:
                nearxy.append(xy[0])
    near_ids[key] = nearxy #gets the coordinates of the snapped point


texts:
I'm trying to create a mask from a feature collection with one feature:
var conus_fc = ee.FeatureCollection('users/jesse/aft-app/conus');
var conus = ee.Image().paint(conus_fc);
Map.addLayer(conus_fc);
Map.addLayer(conus);

The feature collection looks like this:

The painted feature collection looks like this:

I don't know why those gaps are there in the painted image; they are not in the feature collection. I need to get rid of them.


texts: Oil and gas well symbology for geoserver

texts:
I need to symbolise an oil and gas well dataset in geoserver.

I have the dataset symbolised in ArcMap using the standard esri well symbols. Is there an easy way to convert this esri layer file into an SLD? Or, is there some where I can reference well symbology in my SLD using an < ExternalGraphic>< OnlineResource xlink:href.../>??


texts: Calculating compass bearing between two GPS coordinates... with a twist

texts:
I'm building a virtual tour application that uses Google's streetview code to render the panoramas into a non-distorted view.  Each panorama has 1 to n links where the user can select to move on.

I use the GPS coordinates of each location to calculate the compass bearing from the current location to the next.  The issue has to do with StreetView's rendering of the link bearings.

When a panorama loads, one of the data fields required is centerHeading.  This determines the part of the image that is rendered smack dab in the middle of whatever element you're using as a map container.

In addition, whatever value centerHeading has is what StreetView sees as due north.  The link rendering locations are displayed relative to north, which causes quite a bit of inaccuracy with the link bearings.

Check out this for an example. The bearing from the start location to the next is 247 degrees, as calculated from the GPS coordinates of each location.  Due to the way the image was taken, the centerHeading value is 175 degrees, and this throws off the actual links quite a bit. In order to have a correct link bearing, I have to manually replace the calculated bearing of 247 degrees with 0.

Multiply this trial and error bearing testing by over 100 seemingly arbitrary heading definitions, and you just might get a headache.

So my question is, is there a way to factor in the location's heading setting in the link bearing calculations?


texts:
I have a problem with GeoPandas sjoin, I work with Ubuntu and I want to use:
opolskieadresy = geopandas.sjoin(smallopolskie, adresy, how="inner", op='intersects')

and i'm getting error:
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/gis/.local/lib/python3.8/site-packages/geopandas/tools/sjoin.py", line 88, in sjoin
    sindex = right_df.sindex
  File "/home/gis/.local/lib/python3.8/site-packages/geopandas/base.py", line 593, in sindex
    self._generate_sindex()
  File "/home/gis/.local/lib/python3.8/site-packages/geopandas/base.py", line 98, in _generate_sindex
    sindex_cls = get_sindex_class()
  File "/home/gis/.local/lib/python3.8/site-packages/geopandas/sindex.py", line 42, in get_sindex_class
    raise ImportError(
ImportError: Spatial indexes require either `rtree` or `pygeos`. See installation instructions at https://geopandas.org/install.html

When I try to install rtree: sudo pip3 install rtee I'm getting:
Collecting rtree
  Using cached Rtree-0.9.4.tar.gz (62 kB)
    ERROR: Command errored out with exit status 1:
     command: /usr/bin/python3 -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-ib_695_y/rtree/setup.py'"'"'; __file__='"'"'/tmp/pip-install-ib_695_y/rtree/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' egg_info --egg-base /tmp/pip-install-ib_695_y/rtree/pip-egg-info
         cwd: /tmp/pip-install-ib_695_y/rtree/
    Complete output (15 lines):
    Traceback (most recent call last):
      File "<string>", line 1, in <module>
      File "/tmp/pip-install-ib_695_y/rtree/setup.py", line 3, in <module>
        import rtree
      File "/tmp/pip-install-ib_695_y/rtree/rtree/__init__.py", line 1, in <module>
        from .index import Rtree
      File "/tmp/pip-install-ib_695_y/rtree/rtree/index.py", line 6, in <module>
        from . import core
      File "/tmp/pip-install-ib_695_y/rtree/rtree/core.py", line 143, in <module>
        rt.Error_GetLastErrorNum.restype = ctypes.c_int
      File "/usr/lib/python3.8/ctypes/__init__.py", line 386, in __getattr__
        func = self.__getitem__(name)
      File "/usr/lib/python3.8/ctypes/__init__.py", line 391, in __getitem__
        func = self._FuncPtr((name_or_ordinal, self))
    AttributeError: /usr/bin/python3: undefined symbol: Error_GetLastErrorNum
    ----------------------------------------
ERROR: Command errored out with exit status 1: python setup.py egg_info Check the logs for full command output.


texts:
Can someone explain what is the history of BING vs. HERE maps and what are the differences as of today (2015 Q4)?

I'd be especially interested in the Satellite view, as BING maps was one of the providers which allowed OSM referencing. Does HERE maps allow such thing now that they are separate?


texts: Is it possible to export 3D model data from SimUrban?

texts:
Today I came across a 3D modelling package I haven't seen before - SimUrban.

The client wants to migrate their SimUrban project to a more modern system, and the provider who built it is no longer in business.

The person using SimUrban wasn't an expert, and we couldn't see any UI controls to export the data into another format.

I didn't get a chance to inspect the system in detail, but it appeared that there was a directory containing many blah.model files, which presumably comprised the 3D buildings shown in the system.

Does anyone know whether it's possible to export the building models from SimUrban into a commonly-used 3D modelling format such as VRML, DXF, DGB or similar?


texts: ArcGIS for Desktop equivalents to ArcInfo Workstation tolerance types?

texts:
In ArcInfo Workstation (formerly called ARC/INFO), the ARC: command TOLERANCE was used to specify fuzzy, grain, weed, snap and nodesnap as the types of tolerances for coverages.

When creating a feature class, what would be the equivalent of these?

There do not seem to be any options to specify these in ArcGIS for Desktop.