texts:
I'm currently working on MODIS satellite product "MOD16A2" which provides regular 500m's resolution pixel of evapotranspiration. I basically need to extract from the map the mean value of different drainage basins. What I have is the raster map of evapotranspiration with no datum and a sinusoidal projection and shapefile of the basins.

I learnt, from this forum, I should use the zonal statistics choosing as input the raster map and basins shapefile and calculate the average. I'm facing two problems though:

There are many fill values which space from 3276 to 3276.6, which is a problem because they're going to be counted on the average,
I have to reproject the raster into the same datum of the shapefile, but when I try to use the "save as" and select the new CRS it generates another map with different number of columns and rows and also different pixel values.

For the first problem I settled on "save as" as null values all the values in the range 3000-4000 and it apparently works but still I'm going to  have to deal with the 2nd problem.

I also tried to use the raster calculator and use the following formula:
("myraster" < 3276 ) * "myraster" and then using the GDAL conversion to set 0 as NODATA, but now I tried again the raster calculator doesn't transform into 0 the values any more. Also, I have to be able to repeat those steps for about 800 times (the number of maps that I have) to get a time series of evapotranspiration.

I don't have any knowledge on writing scripts yet and I'm only able to use QGIS.

I don't know how to use Google Engine and the time I have left is not much. I understood that reprojecting satellite's products, such as the MOD16 A2 MODIS'S product, can give some problem during the resampling in fact I obtained as result different number of rows and columns and different values. What I did is to use the "save as" right clicking the raster and only set the no data value in the range from 3000 to 4000. In this way the new map, matched perfectly the original one and provides all the values in the range settled before as no data. Later I just changed the projection of the drainage basins into the same projection of my raster, and then applied the zonal statistics. This time seems to work but just with a "small" incovenience: I had to repeat, manually, this step for all the maps. The zonal statistics from GRASS helped me to do a batch process but instead I obtained as many shapefile as the number of rasters is. Since I need the mean value of evapotranspiration which is calculated in table attributes of every shapefile generated, I wonder if there's a way to automatically extract the column I need from every table.


texts: Map service created correctly but is inaccessible from the Flex API

texts:
I created a map service from a mxd document. When I try to display to it through the Flex API:
<esri:ArcGISDynamicMapServiceLayer url="http://arcgisServerIP:8399/arcgis/rest/services/myService/MapServer" />
I get the following response from the server:
{"error":{"code":404,"message":"Service 'myService' of type 'MapServer' does not exist or is inaccessible.","details":[]}}

Does anybody have a solution?

Regards.


texts: Adjusting query in field calculator in QGIS 3.2?

texts:
I have a geopackage file which contains point data with a column of 'depth of cover' and score.

In the score column I am using the field calculation to create an expression so that any number less than 0.5 in the depth of cover column would be represented as a 1 and anything in between 0.6 to 5.0 would be represented as a 2, and anything greater than 5.1 would be a 3 in the score column.

This is the expression I have used:

CASE
  WHEN  "Depth of cover" <=0.5 THEN 1
  WHEN  "Depth of cover" >=0.6 <=5.0 THEN 2
  WHEN  "Depth of cover">=5.1 THEN 3
END

I can get the 1 and 3 values to work but not 2.  I know I have the expression wrong at this point because I'm doing a range rather than a less than or greater.

Following on from this if I wanted to change the 1 to a High 2 to medium and 3 to Low how would I do this as I get Column 'High' not found.


texts: Is there a way to create a new vector layer from a template?

texts:
I would like to be able to set up a template for a vector layer with a standard set of attributes rather the having to specify them time.  Similar to SQL "CREATE TABLE blah LIKE foo"

Is this possible?


texts: Geoserver get UTFGrid in Zoom, X, Y format

texts:
I am trying to overlap an UTF grid on a WMS layer in order to interact with it. As front end library I am using leaflet and the layers are served by GeoServer.

Now, I found some Leaflet plugins that are able to handle the UTFGrids, the problem is that they accept as input an URL formatted in this way:

http:///layer/z/x/y.json

Is there a way to get my UTF grid using this template?


texts:
I was wondering if there was any way to create main axis (axis of symmetry) lines for polygons.
I have polygons that represents buildings.
I will want to find main axes (axis of symmetry) for those buildings and axis direction respect to a cardinal point.
Best regards
Stefano


texts: I have collected the public transport routes and i want to map a map that looks like this. My problem is i am not able to show the shared segments as parallel lines. I have not been able to find good resources on this. I have good experiences with making maps and GIS softwares. If this requires special tools , then i can learn those too. Although i am looking for a static map, any possibility of making an interactive map would be better.


texts: Best programming language for handling both statistical and spatial analysis

texts:
I am frequently involved in research where I do statistical programming and analysis using tabular datasets, and incorporate independent variables that have to be calculated in a GIS.  I currently rely on SAS for statistical programming and the ArcGIS GUI for obtaining the independent variables.  I face 2 challenges with this workflow.  First, I have to export my ArcGIS output and import it into SAS.  Second, I do not have a log of all the procedures I perform in the GIS.  As a result, my workflow is prone to errors because I am transfering data between systems and have no way to look back on my GIS steps.

What are the best programming languages (and environments) for handling both the statistical and GIS analyses in one environment?

Preliminary Thoughts

A partial solution would be to use Python with ArcPy.  I am not a fan of this solution because Python was not designed for statistical programming or statistical analysis.  I am wondering whether I should be looking into R or PostgreSQL with PostGIS.


texts: Error installing QGIS on Mac OS X 10.11.6 (Why do I "need Xcode.app" to "compile qt5"?)

texts:
I'm trying to install QGIS by following the instructions in from this answer (using homebrew).

I get the following error:

$ brew install osgeo/osgeo4mac/qgis3
==> Installing qgis3 from osgeo/osgeo4mac
==> Tapping brewsci/bio
Cloning into '/usr/local/Homebrew/Library/Taps/brewsci/homebrew-bio'...
remote: Enumerating objects: 226, done.
remote: Counting objects: 100% (226/226), done.
remote: Compressing objects: 100% (223/223), done.
remote: Total 226 (delta 1), reused 134 (delta 1), pack-reused 0
Receiving objects: 100% (226/226), 161.66 KiB | 0 bytes/s, done.
Resolving deltas: 100% (1/1), done.
Tapped 211 formulae (245 files, 496KB).
qt5-webkit: Must be built with `brew install --no-sandbox ...`, or install steps will fail.
qt5-webkit: A full installation of Xcode.app is required to compile this software.
Installing just the Command Line Tools is not sufficient.
Xcode can be installed from the App Store.
Error: Unsatisfied requirements failed this build.

Why is it telling me I need Xcode, and why would I need to compile qt5 to load QGIS?


texts: Performing ArcGIS Network Analyst Routing on one side of street before other?

texts:
I am having issues with ArcGIS Network Analyst when I am trying to solve some addresses to find the best route, in a walking perspective.

How can I stop the program from jumping from one side of the street to the other?

I want the street to be walked all odd numbers first, and then all even


texts: ArcPy UpdateCursor appears to not Update field in open Attribute Table?

texts:
I have been trying to get Update Cursor to update a field in my attribute tale.

import arcpy
fields = ['STATE_ABBR','Last_Update']
fc = "C:\Users\sword\Desktop\GIS Web Course\PersonalProject\GIS_Project.gdb\USAStates"
with arcpy.da.UpdateCursor(fc, fields) as cursor:
    for row in cursor:
        if row[0] == 'HI':
            row[1] = '1'
        cursor.updateRow(row)

One Row has state abbreviations and the row Last_Update is blank, I just want to place a one in that row. I enter it in the Python box in ArcGis, it throw no errors but simply does nothing. The row doesn't update. The Last_Update column is a string field, and I have verified the path name in fc. What am I doing wrong? This is almost a copy from ESRI webite of the format.


texts: Map Matching? (Finding Nearest Roads...)

texts:
I am making Pedestrian Navigation System for my school

I created node and arc and digitized buildings with ArcGIS. Then, exported KML to use Javascript Dijkstra's Algorithm to find route.

So far, so good.
The thing is if i want to find the shortest path or fastest route, I have to type NODE number.

I want a person to click (or tab his or her's smartphone screen) a building and automatically snap the point to the nearest nodes (roads)

What keyword should i look up or do you have any idea?


texts:
I do not manage to install QGIS on my windows 10 machine using the osgeo4w installer.

During the installation I get this message:

At the end of the installation I get this message:

And when I subsequently try to run QGIS, I receive this message:

Similar issue can be found in various fora, but none of them provides a conluding answer or did help me to solve the problem.

Any propositions?


texts:
I am a new in this field and need help,  how to code using SQL spatial a function that retrieve loops in a network.The input, flag: Geography and barrier: Geography  and   two spatial tables, first one will list the edges(edgekey, from_junction, to_junction, polyline) and the second table containing juctions(junctionkey, point).

The output will be a spatial table containing (edgekey, polykey,  from_junction, to_junction, polyline), a single polygon  will have the same polykey for each from_junction  to_junction value.


texts: PostGIS retrieve polygons within a rectangle

texts:
I am building a system to bring in shapefile data from a PostGIS database and render it on top of OpenStreetMap tiles using the Leaflet.js library.

My table (sierra_county_parcels) has a polygon geometry field (geom) in SRID 4326, and a field called "upc" that contains the parcel ID. There are about 16,000 rows in sierra_county_parcels. What I need to do is retrieve all of the upc values for rows where the geom polygon is within the currently visible part of the Leaflet map control. I have already written the code to send the northwest and southeast lat/lon pairs and generate an SQL query, which looks like this:

SELECT upc AS pt_parcel_key
FROM  sierra_county_parcels
WHERE  ST_Within(geom, ST_Envelope(ST_GeomFromText('MULTIPOINT(33.167444534375925 -107.14931488037108, 33.08909661229697 -107.35513687133788)')))

The MULTIPOINT WKT is what I'm currently using to specify the corners (nw_latitude, nw_longitude, se_latitude, se_longitude).

The problem is that the query always returns 0 results. For this query, the map was zoomed out to show most of the region where the geometries should be.

Any ideas?

Thanks in advance,

John W.


texts:
I'm trying to create an indoor map(from floor plan) in GeoServer and display in via OpenLayers.

The only way that I found is to create a Store in GeoServer(in my case, using ImagePyramid) from georeferenced .tiff, which one has a local bounding box([0, 0, 11, 9] in EPSG:3857, in my case).

It works fine, but is there another way?

I'm looking to custom CRS, but I don't know whether it's worth the effort.


texts:
I am trying to find the shortest path between two points based on an attribute other than length (like speed limit). I am using ArcMap 10.2.1 (no network analyst extension) and FME 2014; but I am trying to just use FME for my analysis.

Scenario:
Each road segment has a length and speed limit, need to find shortest path from A to B while travelling under a certain speed.

Data: road network

So right now I am filtering out the road segments in ArcMap based on the speed limit (i.e. speed limit < 60km), and exporting that out as a shapefile, then using that as the input network. In FME I am using CREATOR to make the start/end points and SHORTESTPATHFINDER to calculate the shortest path based on the filtered network and length of the road segments. Are there any ways to just start with the original dataset and set my criteria for attributes (speed limit and lengths of segments) all in FME?


texts:
I have installed the version 2.18.24 in QGIS and apparently the plugin 'Digitizing tools' does not include the feture 'add ring' but it was replaced by 'fill ring with a new feature'. However, somehow I cannot create this new feature neither.
Case this would work... if I want a hole in the polygon, do I have to create a feature and then delete it?
Do anyone have some information on that?


texts:
I have wms tile layer that only show street name labels. An ajax function pulls in the street line json layer to the map. The line get drawn on top of them labels. Is there a way of controlling z-index of the json layer to be drawn under the tile layer?

I have tried .addTo(map).setZIndex(9); on my wms layer but the json layer still get drawn on top.


texts: How to calculate whether a point is in a shape using geodetic calculations

texts:
I’m new to GIS so forgive me if my terminology isn’t completely correct.  Really liking the options available in the GeoTools library but I’ve hit a stumbling block that I’m hoping someone will be able to advise on.  I’ve looked at the classes in the library but can’t find a method that is giving me what I need.  I have a scenario where I have a point and a polygon defined by latitude and longitude (WGS84).  I need to determine whether the point is in (or on the boundary) of the polygon.

I have tried using the Intersects method but it looks to be treating the polygon as if it is on a 2D plane.  The bottom left and right corners of my polygon have the same latitude value but they are spaced longitudinally by 10 degrees, so it’s quite a long edge.  There are no intermediate points along this edge of the polygon.  The Intersects method appears to use a straight line between these two points, effectively taking this bottom edge under ground when on a globe.

What I need is for this polygon’s edges to be interpreted as great circle lines between points on the surface of the globe, thereby creating slight curves, so when we determine whether a point is in the polygon, we are using the shape with slightly curved edges.  Obviously the greater distance between polygon points the greater the curve.

Is there a method to achieve this in the GeoTools library?  If not, does anyone have alternative open source Java libraries that can solve what must be a common request?


texts: Using ArcGIS Desktop on a remote desktop server

texts:
We have about 5 users, they will user ArcGIS Desktop 10.2 by connecting to a Server through Remote Desktop Connection. ArcGIS Desktop obtain licence from a licence manager on anther server.
If the 5 users working at the same (Remotely from same Server), how many licences on Licence manager will be used?
Are the licences are per IP or per user?
Adding one more point
How I can give the users different type of licences? Example one will have Advanced and other will have Editor, etc..
We tried to experiment this, If licence changed for one user it will change for all other user.


texts:
I'm building an application that uses data from GPS units on board vehicles and displays it on a map. A lat/long is recorded every .5 seconds and a session can be many hours... this can add up to be a lot of data. The vehicles travel in a straight line the majority of the time but the data needs to be very accurate when a vehicle turns a corner, hence the high recording frequency. I'm storing the lat/longs as a POINT in a PostGIS/PostgreSQL database (This isn't set in stone).

I'm trying to find a way of reducing the dataset when displaying it on the map while still maintaining accuracy.

I have come up with a few options:

Group the data into 2 second blocks and average the lat/longs in each group. This is really easy to do but has too much information when traveling in straight lines and not enough when turning corners.
Group points together which are a certain distance (a vehicles length?) and time (4 seconds?... to ensure it won't be group if a vehicle moves over the same point half an hour later)
Calculate splines of the data (I don't know very much about this option).

Perhaps PostGIS already provides this functionality? I'm simply looking for someone to point me in the right direction.


texts: Using shortestpathpointtopoint tool with PyQGIS

texts:
I am a total beginner with Python and programming in GIS in general. I have point data and I want to compare in which month people walked the longest distance.
It is a relative large dataset and I want to do it relatively automatically.
This is my approach:
import processing

layer_anfang = QgsProject.instance().mapLayersByName("MVG_2016_Wetter_Januar")[0]
layer_ende = QgsProject.instance().mapLayersByName("MVG_2016_wetter_januar_ende")[0]

for i in range(0, 5):
    startlat = layer_anfang.getFeature(i).attribute(3)
    startlon = layer_anfang.getFeature(i).attribute(4)
    endlat = layer_ende.getFeature(i).attribute(5)
    endlon = layer_ende.getFeature(i).attribute(6)

    x = startlon, startlat, endlon, endlat
    print(x)

processing.run("native:shortestpathpointtopoint", \
{'INPUT':'C:\\Users\\.......', \
'STRATEGY':0,'DIRECTION_FIELD':'','VALUE_FORWARD':'','VALUE_BACKWARD':'','VALUE_BOTH':'','DEFAULT_DIRECTION':2,'SPEED_FIELD':'','DEFAULT_SPEED':50,'TOLERANCE':0, \
'START_POINT':'startlon, startlat [EPSG:4326]','END_POINT':' 11.47498, 48.10944 [EPSG:4326]',\
'OUTPUT':'C:/Users/.....'})

The plan is to iterate over the attribute table and get the lon/lat data from my attribute table and then parse it into the parameters for "START POINT" and "END_POINT". As you can see in this example I did it for the Start_Point, the END_Point is as it was before. I always get the error wrong parameter for START_Point, even though, if I print it for example it is in the right style.
Must I cast any datatype or something on the variable or is this all complete wrong?


texts: Pros and cons of saving data on a memory stick

texts:
I have just about reached the capacity of my laptops internal Solid State Hard Drive (Sandisk SD8SNAT-256G-1006) and I looking for some additional storage without having to buy a new machine.

You can now buy flash drive's with 256 gb of memory for next to nothing. My thoughts are to store all my data including shapefiles, PostGIS and SpatialLite databases on a flash drive which will also be useful if I want to work with the data on another machine.

Can people please advise on the downsides to this kind of set up? Will QGIS be able to process the data of a cheap external at the same rate/ efficiency as it can on my existing internal?


texts: Classifying type of land use each spectral represents from remote sensed data?

texts:
How do you know the type of land use each spectral represents from an remote sensed data into a classification ?

Which classification is the most typical to use when doing that ?

How many classification are there do you recommend ?


texts: Convert Longitude/Latitude to GTA V coordinates

texts:
This problem is very similar to this post from awhile back, but it was never solved. I am trying to create an interactive GTA V Map with the Leaflet JS API. I want to convert the longitude and latitude values on the map to GTA V coordinates. Someone replied to other post saying I would need to make the map use a cartesian coordinate system and linked this page from the Leaflet documentation. I still do not understand how to apply this to my question. Here is a coordinate map of the GTA V map.

Can anyone walk me through this or point me in the right direction?


texts: Is there any faster ways to count intersections between lines and polygons for large datasets with python?

texts: I utilized the Polygon.inetersects() method to judge intersections between line dataset and polygon dataset with "shapely", but it turn out to be slow but I have a large number of data. Can anybody give any suggestions?


texts: How to add Leaflet controls in an outer div

texts:
What is wrong with this code ?

var mapbuttons_div = L.DomUtil.create('div', 'mapbuttons leaflet-bar leaflet-control');
var zoom = L.control.zoom(
    {
        position: 'topleft'
    }
);

mapbuttons_div.appendChild(zoom.getContainer()); // error

I want to add the zoom buttons div to an outer div called mapbuttons. But I get an error:

  Uncaught TypeError: Failed to execute 'appendChild' on 'Node': parameter 1 is not of type 'Node'.


texts: What is the problem since i have "test_script" folder  extant in C: but arcpy detects as a truncated one (i.e est_script) and signals error.


texts: Using OpenLayers 3 forEachFeatureAtPixel to specify layer?

texts:
I must first state that my skills at javascript are about the same as my skills in Spanish. I know barely enough to get myself into trouble and then I go looking for help. anyway I have an application that has a few geoJSON layers in it. I have created a map.forEachFeatureAtPixel function to update a div with some information.

The problem is it looks at all of the layers and I just want it to look at one. this way I can create another function to get some different attributes from another layer. I have looked over the documentation and I see the opt_layerFilter option but I dont know how to apply it (due to my lack of Javascript knowledge). If someone could give me a hint here hopefully I will then be able to read the documentation and stop asking these questions.

function:

var displayFeatureInfo = function(pixel ) {
  var feature = map.forEachFeatureAtPixel(pixel, function(feature, layer) {
    return feature;
    });

  var info = document.getElementById('info');
  if (feature) {
    info.innerHTML =  '

' + feature.get('Name') +
     '

' + '
API:

' + feature.get('API') + '
'
     + '
Status:

' + feature.get('Status') +  '
' +
     '
Legal :

' + feature.get('Legal') +  '
' +
     '
Site Name:

' + feature.get('Site_Name') +  '
'+'
';
  } else {
    info.innerHTML = '&nbsp;';
  }

  if (feature !== highlight) {
    if (highlight) {
      featureOverlay.removeFeature(highlight);
    }
    if (feature) {
      featureOverlay.addFeature(feature);
    }
    highlight = feature;
  }

};


texts: How to integrate Proj4 into an iOS / iPhone project?

texts:
This may be overkill, but I need to reproject a point coming from the CLLocationManager, which is in EPSG:4326, to Dutch RD (Double Stereographic, EPSG:28992) in my iOS application and I am using Proj4 to do it. At least that's what I set out to do. Here's what I tried:

I grabbed a Proj4 Xcode project from route-me and added it to my Xcode project.
I added a dependency on Proj4 to my build target. (Target -> Properties -> General -> Direct Dependencies -> + (add)
I added the path to the Proj4 project directory to my Library search paths. (Target -> Properties -> Build -> Library Search Paths)
I added the path to the Proj4 project directory to my Header search paths. (Target -> Properties -> Build -> Header Search Paths)
I added a reference to Proj4 (-lProj4) to my linker flags (Target -> Properties -> Build -> Other Linker Flags)

But I still get a linker error when compiling:

ld: library not found for -lProj4

What am I missing here?


texts: Critical facility pipe connections in python

texts:
I'm trying to find a way to determine pipeline connections located downstream a critical facility and score how far those connections are from the facility.

Is there a way to accomplish the above task with python? I have a shapefile and would like to work with that to write a script in python.


texts: Setting coordinate Grid with map composer in QGIS

texts:
I have been trying to display a coordinate grid on my map with QGIS.
I have a shapefile projected to UTM zone 34N. I have set the CRS to "WGS 84 Zone 34N" on the main properties tab and the same for map grid properties on the CRS tab. My interval is map units. Unfortunately, I am definitely missing or doing something wrong and I get some weird lines at the corners (see image) of my map and the same when I change the interval to cm.

My main goal is to display normal coordinate grid on my map. If you can guide me on this or any input about the process of doing this correctly.


texts:
The EPA has an implementation called How's my water. https://mywaterway.epa.gov/
As I understand it, here is kind of how their application works.

Input -> lat/long or address in the USA
The application gets the 12 digit HUC associated with the location e.g. "080801030109"
Using the HUC and a webservice called ATTAINS the application gets the assessment Unit Id e.g. "LA060801_00"
Using the assessment Unit Id and ATTAINS, the application returns useful data e.g. "the name of the nearest waterbody e.g. "Vermilion River-From headwaters to LA-3073 bridge"

ATTAINS -> webservice ATTAINS
So my question is how do I approach step 2 above?
Can Google Earth Engine return the 12 digit HUC associated with a lat/long or address?
here's the code so far and the result.


texts: How do I add multiple indices as bands to an image?

texts:
I want to be able to classify my image using bands for prediction (B2, B3, B4, etc.) but I also want to use indices as "bands" (B2, B3, B4, NDVI, EBBI, MSAVI2, etc.). Unfortunately, I can't seem to think of the proper script for this and it keeps giving me the error:

Here is my script:

// Creating a buffer around my region of interest
var buffer = ee.Geometry.Point([103.83841283715412,1.3736828122619071]).buffer(10000);
Map.centerObject(buffer, 12);

// Load Landsat 8 surface reflectance data
// Look for adequate image filtered by date and cloud cover and region
var l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
              .filterBounds(buffer)
              .filter(ee.Filter.calendarRange(2018,2019,'year'))
              .filter(ee.Filter.calendarRange(4,5,'month'))
              .map(function(image){return image.clip(buffer)});

// Function to cloud mask from the Fmask band of Landsat 8 SR data.
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = ee.Number(2).pow(3).int();
  var cloudsBitMask = ee.Number(2).pow(5).int();

  // Get the pixel QA band.
  var qa = image.select('pixel_qa');

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

  // Return the masked image, scaled to [0, 1].
  return image.updateMask(mask).divide(10000);
}

// Map the function over years of data and take the median.
var composite = l8sr.map(maskL8sr)
                    .median();

//The Band Variables for Landsat 8
var B2_Blue = composite.select('B2');
var B3_Green = composite.select('B3');
var B4_Red = composite.select('B4');
var B5_NIR = composite.select('B5');
var B6_SWIR1 = composite.select('B6');
var B7_SWIR2 = composite.select('B7');
var B10_TIR1 = composite.select('B10');
var B11_TIR2 = composite.select('B11');

// Creating an add variable function for Landsat 8 NDVI.
var addNDVI_L8 = function(image) {
  var ndvi = composite.expression(
  '(B5-B4)/(B5+B4)', {
    'B5': B5_NIR,
    'B4': B4_Red
  }).rename('ndvi');
  return composite.addBands(ndvi);
};

// Creating an add variable function for Landsat 8 MSAVI2.
var addMSAVI2_L8 = function(image) {
  var msavi2 = composite.expression(
  '(2 * B5 + 1 - sqrt(pow((2 * B5 + 1), 2) - 8 * (B5 - B4)) ) / 2', {
    'B5': B5_NIR,
    'B4': B4_Red
  }).rename('msavi2');
  return composite.addBands(msavi2);
};

// Creating an add variable function for Landsat 8 EBBI.
var addEBBI_L8 = function(image) {
  var ebbi = composite.expression(
    '(B6-B5)/(sqrt(B6+B10)*(10))', {
    'B5': B5_NIR,
    'B6': B6_SWIR1,
    'B10': B10_TIR1
  }).rename('ebbi');
  return composite.addBands(ebbi);
};

//Defining the composite image with the added NDVI band
var composite_NDVI = addNDVI_L8(composite);

//printing to see if it has 13 bands (i.e. with the NDVI band added) (THIS WORKS!)
print(composite_NDVI, 'Image with added NDVI band')

//THIS CODE DOESN'T WORK "composite.map is not a function"
var composite_indices = composite.map(addNDVI_L8).map(addMSAVI2_L8).map(addEBBI_L8);

print(composite_indices, 'composite indices');

so my script, doesn't seem to work:

var composite_indices = composite.map(addNDVI_L8).map(addMSAVI2_L8).map(addEBBI_L8);


texts: Google Earth Engine_Conversion to vector

texts:
I have a classified geotiff raster that has two classes 0 and 1. I changed its type to int and imported it from my computer to GEE using upload a new asset and I changed its var name to image and I wrote this lines:
Map.addLayer(image);

var classes = image.reduceToVectors({
      reducer: ee.Reducer.countEvery(),
      scale: 20,
      maxPixels: 1e10
});
var result = ee.FeatureCollection(classes);
Map.addLayer(result);

// Export the results to Drive
Export.image.toDrive({
      image: result,
      description: 'SGenova_int',
      scale: 20,
      crs: 'EPSG:4326',
      maxPixels: 1e10
});

I got this error while exporting:

Error: Image.setDefaultProjection, argument 'image': Invalid type. Expected type: Image. Actual type: FeatureCollection.


texts: IRasterLayerExport3 doesn't save raster until ArcMap closes

texts: I'm trying to export raster using IRasterLayerExport3, but it saves raster file only after closing ArcMap. It seams that IRasterLayerExport3 creates placeholders for raster (empty file - 3kb), and it keeps all information in memory or maybe in temporary file. And it writes data to an empty file only after creating a new map document or closing ArcMap. Any ideas about this?


texts: Select overlapping features from single layer

texts:
in QGIS 3.10 LTR we have the opportunity to identify vector features with the mode "Top down".
Is there a way to do this with the "Select features"? We'd like to select all overlapping features when using the single click.


texts: Incorrect field name resulted in RuntimeError: ERROR 999999 (Also, better code for adding and calculating fields)

texts:
This is the bottom half of my code (and it's been edited somewhat to share publicly). The code works perfectly fine when the fieldList has only three fields in it, but when I populate it with all 108 fields, I receive the following error, just before it starts to calculate fields:

Traceback (most recent call last):
  File "C:\GIS\1GIS_DATA\Folder1\Folder2\Python and Tools\RunThisPythonCode.py", line 125, in <module>
    value = row.getValue(field)* row.anotherField #<--This is the line the error is referring to
  File "C:\Program Files\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects\arcobjects.py", line 944, in getValue
    return convertArcObjectToPythonObject(self._arc_object.GetValue(*gp_fixargs(args)))
RuntimeError: ERROR 999999: Error executing function.

Here's the code:

....

fc = arcpy.CopyFeatures_management(layerName, outFeatures)

#---------------------------------------------
# Add fields to a feature class:
#---------------------------------------------
fieldList = ["field1", "field2", "field3", ............., "field99"]  #<-- fieldList
fieldListx = []
for f in fieldList:
    fieldx = f+"x"
    arcpy.AddField_management(fc, fieldx, "DOUBLE")
    fieldListx.append(fieldx)

#---------------------------------------------
# Calculate fields:
#---------------------------------------------
rows = arcpy.UpdateCursor(fc)
for row in rows:
    for fieldx in fieldListx:
        field = fieldx.rstrip('x')
        value = row.getValue(field)* row.anotherField  #<--This is the line the error is referring to
        row.setValue(fieldx,value)
    rows.updateRow(row)
del row, rows

Is this an issue of memory or something? Why does it produce an error before Calculate Fields when I add the full list of 108 fields, but not when I run a small sample of only 3? Thanks!!


texts: Find out point which falls between two parallel lines

texts:
I am facing one problem in ArcGIS. I work on a navigation database. In our database single lane streets are represented by a single Line, while a multi-lane street (street with divider in centre) is represented by two parallel lines (red coloured lines in picture).

I have a point shapefile with some points falling inside multi-Lane street and some outside.

I want to create an ArcPy script which would find the points which fall inside Multi-Lane Streets. i.e. between these parallel lines (marked in picture).

I do not know how to achieve this, Can somebody help me?

I did some exercise on it and i found that creating buffer on one side of line can create inside Multi-Lane polygon(Shown in Picture).

but now the problem is , polygon is actually crossing the line(i.e. overlapping the multi-lane boundary). so it will catch unnecessary points. is there any way to align this polygon to the street line ?

Note: integrate will not work here, because it also moves the street line. i need to just align polygon along to the street line.


texts: CSV matrix of coordinates to Raster in QGIS

texts:
Is it possible to create a raster in QGIS from a CSV file that is a “matrix of coordinates” without transforming it into a vector?

I explain: My data is a matrix of 5500 columns and 175 rows, it is 962500 points. The Columns are the X coordinates and the Rows are the Y coordinates. The values in the matrix are the Z coordinate of each point. The X distance between points is 0,04 cm and the Y distance is 0,80 cm. I could try reorganizing the data in a way that I have just three coordinates in columns X, Y, Z for all the 962500 points. However, it means an additional step and I would have to do it repeatedly as my experiment produces the CSV file described before.

How can I create the raster in QGIS from the original matrix?


texts: Finding "most frequently intersected area"?

texts:
I have a geodataframe that contains geometries of droughted area for several years. To identify consistently droughted areas, I want to find areas that are droughted for the largest number of years.
In the image below, the 3 shapes represent droughted area outlines over 3 separate years. I want to end up with a geodataframe (or list of shapely polygons) with a count of the number of overlaps. So area A = 2, B = 3, and C = 2. I can then symbolize by color to show areas that are consistently droughted.

I've tried to get intersected areas only then perform a spatial join, but I can't figure out how to perform a "count" using GeoPandas. This answer seems close to what I want, but I wasn't able to translate it to Python (not using ArcPy).


texts:
My current goal is to calculate the following
(a<10)*a and (b<0.7)*b and c=2 and d=2
this formula already worked in the native raster-calculator and now i'm trying to get this to work within the modeler. Long story short i can't get it to work.

got it to work with ((a<10)*a+(b<0.7)*b+(c=2)+(d=2)) but it's not the result i was hoping for, with the native calculator i get a new raster where one Attribute is NULL and the other one has a value which represents the area i'm looking for and now i get too much information...

is there anyway to include logical operators like 'and' to achieve this?


texts: Generating ESRI compatible raster reference files in QGIS?

texts:
Backstory: employer is GIS incapable and most of the data we work with doesnt play nice with GIS software. That said, we have a client who is very GIS oriented. My boss (not my work) is providing them with what are really just pictures of an environmental variable that he creates in photoshop. They are asking for the "world file" for the rasters so they can plug the images into their ESRI powered dataflow.

I have georeferenced the rasters before to make some overlays, but I believe they are currently referenced to NAD (it's the Gulf of Mexico-NAD worked).

I have seen some guide info online for generating world files, projection files, etc with QGIS but wanted to ask the community's opinion on how to make the simplest, most compatible solution? w

ould it be the '.wld' extension I've seen or do i need the projection '.prj' file?

Do I need to find out the clients projection since they are likely on a more global coordinate system than the NAD?


texts: calculate total Miles in a polygon using Google Map

texts: I have an application that allows the users to create any custom polygon layer on top of google maps, i am interested to know if their is a way to select the road layer in this polygon and then calculate the total miles for these roads.


texts: unable to create or edit any layer on geoserver

texts:
I've been using Geoserver for some time now with no problems until recently.

After adding some layers and working on my layouts now I am unable to add or modify layers.

I can add stores and create workspaces but if I try to add a new layer, after clicking the save button the server just reloads in the same "Edit Layer" page and no changes are saved nor a layer is added.

I've tried to reboot, uninstall, install windows service, WAR, different version and got the same result. Then I tried to restore my Geoserver data dir from my backup and layers are back but I can not modify nor add new ones.

Geoserver logs show nothing unusual, no errors, no warnings or anything else to provide any details.

May it be some windows configuration that went bad?
I´m on windows server 2008


texts: Pasting more than one geometry at a time in Quick WKT plugin for QGIS

texts:
I have a CSV of WKBs that I need to convert to WKT.

Can I paste them all into Quick WKT at the same time?

How do you break the spaces between geometries?

Page break, comma, semicolon, etc. do not seem to work.


texts: Calculating mean of a field using Python

texts:
I am very new to ArcGIS and Python. I want to do a basic calculation and find the mean of 1 field of a table.
The table is called Top 10 Tracks and the field in it called Path Length (km)
and there are 10 values in that field and I want to find mean of those values.
I have been googling but I couldn't find any information about this (maybe I wasn't googling enough)
Could you please help me write the python code for that or recommend a material I can read?


texts:
I recently been using FME to publish feature layers but have now had to move to map image layers because of the issue around grouping limitations with feature layers.
I am now trying to work out how is the best way once I update a geodatabase or a enterprise database with data, how it can be automatically published as a map service? Is there a straight forward process to do this?


texts: PyQGIS export to JSON: writeAsVectorFormat and precision

texts:
In my QGIS plugin qgis2leaf I am exporting all layers as JSON into a JSON directory. Therefore I would like to enable the user to define the precision in the GUI. The normal "save as" dialogue in QGIS has the possibility to define the precision:

Is there any possibility to use a precision in the writeAsVectorFormat?
The found documentation gives me no clue:

QgsVectorFileWriter::WriterError QgsVectorFileWriter::writeAsVectorFormat   (   QgsVectorLayer *    layer,
        const QString &     fileName,
        const QString &     fileEncoding,
        const QgsCoordinateReferenceSystem *    destCRS,
        const QString &     driverName = "ESRI Shapefile",
        bool    onlySelected = false,
        QString *   errorMessage = 0,
        const QStringList &     datasourceOptions = QStringList(),
        const QStringList &     layerOptions = QStringList(),
        bool    skipAttributeCreation = false,
        QString *   newFilename = 0
    )       [static]


texts: Looking for compute intensive GIS vector operations

texts: I am looking to make efficient algorithms for those vector GIS operation that take a long time for GIS software like Arc-GIS, Grass etc. I have considered polygon overlay but looking for other operations as well.


texts: I have a problem, which is, I have a table ways (like in the Workshop of pgRouting) and I also have a different table street_name, where it contains the name of the street (among other information). Both have a column geom and I want to connect a given edge of the pgr_dijkstra calculated route with the name of the nearest point in street_name. I'm already using a version of this function to convert lat/long in source/target values. Is it possible?


texts: gt-geojson not able to read featurecollection


texts: Creating sidebar help information for Python script tools of ArcGIS Desktop?

texts:
I have been creating custom tools (Python script tools) using Python with the ArcPy module.

I can set up the tool successfully with parameters and all, but I cannot figure out how to create help text that shows up on the sidebar of the tool like the standard ArcGIS tools have?


texts: Converting raster layer to simple feature in R

texts:
I have read my raster data using

r <- raster("F:\\Studies\\Whole2015.tif")
class(r)
> class(r)
[1] "RasterLayer"
attr(,"package")
[1] "raster"

How can I transform my data to simple feature?

I want to use a function named st_geometry, which can only be used to sf object.


texts: How can I generate a specific number of random points in a polygon with PostGIS?

texts:
I have a postgis table that I've uploaded a shapefile to, this shapefile contains 200+ polygon rows with several integer columns.

What I'd like to do is similar to the QGIS 'random points in polygon' function except with numbers that are much larger. My initial approach with QGIS crashed upon trying to add ~500,000 random points to each polygon - I need the number of random points to come from the related integer column.

After some research, an indexed postgis database is the way to go. But I'm unsure how to generate random points in a polygon for each polygon in my table, and to have the number of points that get randomly generated, come from a corresponding column.

New to Postgres, and postgis, but familiar with mysql, QGIS, ArcGIS, and geographic concepts.

EDIT:

Here is a picture to better show what I would like to accomplish. Each polygon has an associated PNTSneeded value, I've labelled the polygons in the screenshot with that number (also shown in the attribute table). I need to generate that amount of points, randomly dispersed throughout each polygon row.

I would only need to do this once, so it wouldn't be an on-going process, although I may need to do the same process for different sets of polygons in the future (I'm sure I could recreate the steps).

In the event that this is indeed too many points to generate, I could always run it twice at half the numbers and then just merge the files into one (might add an extra level of randomness). But I would much rather be able to do this in one go.


texts: Starting Dash Space Pattern with space instead of dash in QGIS?

texts: I need to create a simple line symbol which starts with a gap (space) instead of a dash. I avoid replacing the simple line with a marker line because of the x-axis offset that i need to introduce. Can be created an expression to realise this type of dash pattern?


texts: Changing displayed rasters in GeoServer?

texts: I loaded to CoverageStore -> ImageMosaic 3 rasters for different date. On viewing layers I need to be able select what tif I viewing now, but I do not see way to to select what I am viewing now.


texts: relatedQuery.setDefinitionExpression not working jsapi arcgis

texts:
I want to relate two tables STATION and TEMP this my code :

var relatedQuery = new RelationshipQuery();
relatedQuery.objectIds = [graphic.attributes.objectid];
relatedQuery.relationshipId = 0;
relatedQuery.outFields = ["*"];
relatedQuery.setDefinitionExpression = "CAST(date_ AS CHAR(8)) LIKE '201201%'";
console.log(relatedQuery);[![enter image description here][1]][1]
station_lyr.queryRelatedFeatures(relatedQuery, function(relatedRecords){
for(var oid in relatedRecords) {
console.log("related :" + relatedRecords[oid]);
}}

my problem is that the properties setDefinitionExpression not working i have always features: Array[366] ( records : 366 )

but when i used this parametres in REST API INTERFACE i find records : 31


texts: In my QGIS-version (2.18.12, Windows) when opening the expression string builder or the field calculator there is no inline help for the functions. I actually never managed to activate this. Is there an option that I am missing?


texts: GEE Error: Provide 'geometry' parameter when aggregating over an unbounded image

texts:
I'm trying to calculate the mean precipitation for a given polygon for a specific period of time. My end goal is to get one polygon with one value that will be the mean precipitation for a given polygon in specific area. The idea is to create on the end table with different years and on each year to be able to say how much rainfall it had between October to march.

the logic of the code, which may be wrong, was to acess the dataset from PERSIAN, to clip it according to my polygon, then to sum all the images values from the months 10-03, and then to calculate mean using ReduceRegion.

It may be wrong logic so if you have any better ideas let me know.

I get error when I try to calculate the mean:
"Dictionary (Error)
Image.reduceRegion: Provide 'geometry' parameter when aggregating over an unbounded image."

why am I getting this messege?

this is my code:

var dataset = ee.ImageCollection('NOAA/PERSIANN-CDR')
                  .filter(ee.Filter.date('1999-10-01', '2000-03-31'));
var precipitation = dataset.select('precipitation');

var clippedCol=dataset.map(function(im){
  return im.clip(geometry);
});

// //First try to get the mean precipitation for a given area

var summing=clippedCol.sum();

// Reduce the region. The region parameter is the Feature geometry.
var meanDictionary =summing.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry: summing.geometry(),
  scale: 30,
  maxPixels: 1e9
});

// Get the number of images.
var count = dataset.size();
print('Count: ',count);
print(clippedCol);
print(dataset,'dataset');

////create list and check some images
// var listOfImages =(clippedCol.toList(clippedCol.size()));
// var firstImage = ee.Image(listOfImages.get(99));

var precipitationVis = {
  min: 0.0,
  max: 200.0,
  palette: ['ff1105','fbff09','28ff25','03fff3','3907ff'],
};

Map.centerObject(geometry);
Map.addLayer(summing, precipitationVis, 'sum');
print(meanDictionary);


texts: Regular line pattern of flat pixels after using gdalwarp

texts:
I observe a regular pattern of flat pixel bands on a DEM reprojected using gdalwarp:

The red line represents a profile for which the following graph is an extract between the two squared dots, so it crosses two flat pixel bands that are clearly visible:

I simply reprojected an original WGS84 DEM to the UTM of the region with:
$ gdalwarp -t_srs EPSG:23032 input_dem.tif reprojected.tif;

The green line is the profile of the original WGS84 DEM, it's super smooth and it doesn't show the flat pixels pattern.
Is this a known issue? And what could I do to avoid this pattern?
Info:
The version of GDAL I use is from this Docker image ubuntu-small-latest available at https://hub.docker.com/r/osgeo/gdal
# gdalinfo --version
GDAL 3.3.0dev-33cf0e31a992be112b3091f012368d15605ed51d, released 2021/03/17


texts:
I have a geojson file with the location of some universities. This is the map.
The logos are stored in a folder named img Here are the logos.
Now I want to replace the default markers with logo for each point.

I could do it one by one but that's not what I'm looking for.


texts:
I am trying to execute standalone python script outside QGIS. But when I try to initialize Qgs Application , it give me below error . I don't have any idea how to get rid of it .

from qgis.core import *
QgsApplication.setPrefixPath("C:\OSGeo4W64\apps\qgis")
qgs = QgsApplication([], True)

  qt.qpa.plugin: Could not find the Qt platform plugin "windows" in ""
  This application failed to start because no Qt platform plugin could
  be    initialized. Reinstalling the application may fix this problem.


texts: How to attribute polygon with names of polylines it contains?

texts:
I have a set of polyline shapefiles and a polygon shapefile. The task is to split the polylines by the polygon border and write the names of polyline objects located within the polygon to the attribute table of the polygon shapefile. It's clear to me how to split lines by a polygon; the problem is how to get the names and add them to the attributes of the polygon shapefile.

Is there any way to do it in ArcGIS (ArcView license)?


texts: ST_AsMVT single layer with multiple Id's

texts:
I'm trying to create mvt using postgis by following this tutorial

If understand correctly ST_AsMVT is used to create a layer and generated layer will have id. same id can be used to apply styles.

SELECT ST_AsMVT(q, 'buildings', 4096, 'geom')
FROM (
  SELECT id, name, type,
    ST_AsMvtGeom(
      geometry,
      BBox(16597, 11273, 15),
      4096,
      256,
      true
    ) AS geom
  FROM import.osm_buildings
  WHERE geometry && BBox(16597, 11273, 15)
  AND ST_Intersects(geometry, BBox(16597, 11273, 15))
) AS q;

Will create a layer named admin with id admin. styling can be applied by using this id - from index.html

{
  "id": "building",
  "source": "osm",
  "source-layer": "buildings",
  "type": "fill",
  "paint": {
    "fill-color": "#D8CEC3",
    "fill-outline-color": "#C5C0BB",
    "fill-opacity": 1
  }
},

How can we create a single layer with multiple id's so we can style them separately.

I'm complete in newbie to PostGIS and Tiles stuff any other approaches and suggestions are welcome


texts: Finding variable names for SF1 Census files?


texts: Will a tile cache serving PNG tiles be slower than JPEG?

texts: I'm trying to figure out the best format to serve my tiles. I'm creating an ESRI tile cache. It will be consumed internally through ArcMap from our server, not online. There are a lot of different options. From what I've read, PNG tiles are much larger than JPEG, but have better quality. Does that larger size produce a noticeable different when serving the PNG vs JPEG tiles, or will there be no difference? I'm looking for good quality and performance.


texts:
So continuing with my last question (Create 3D dxf with OGR) I have decided to give the GRASS option a go.

I have imported my polylines into GRASS using v.in.ogr and converted them to 3D using v.to.3d.  I checked the lines using the GRASS info tool and the polylines all have the correct z values.   I then exported it using v.out.dxf but this is where I am having troubles.

The dxf contains the z values for the polylines:

Extract

VERTEX
8
line
10
401120.708496
20
6875185.007069
 30
539.000000
0

However when I open it in AutoCAD Civil 3D it doesn't recognize the z value and just makes all the lines 2D.  I have checked them with the info tool in AutoCAD but the z values are all the same.

Is the dxf created using v.out.dxf supported in AutoCAD in order to get correct 3D?


texts: "Export map to image" with a transparent background in QGIS

texts:
I need to export (a lot of) images from QGIS using the "export map to image" command and I need them to have a transparent background.

Is that possible, and how could I achieve it?


texts: I am using ArcMap 10.3. In my map there are resources, and I want to find the location which minimize the distance between them. If I use the Mean or Median Center weighted on the resource I find just one point. Is there any tool in which I can find 2, 3 points which minimize the distance?


texts: I'm trying to find a way to determine pipeline connections located downstream a critical facility and score how far those connections are from the facility. Any Thoughts?


texts: How to symbolize WWF Ecoregions in different colors in QGIS?

texts: Why is a shapefile shown as coloured independent features in ArcGIS shown in a single colour in Quantum GIS? What can I do to display this layer (wwf ecoregions) with the different ecoregions specially coloured?


texts: Trying ST_Slope with SRTM rasters and a shapefile but the output raster has different SRID

texts:
I have two SRTM rasters https://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_06/ (N40W074 and N40W075) with Datum "WGS_1984"

I've united them to be able to work with one raster:

CREATE TABLE rasters.srtm AS
SELECT ST_Union(rast, 1) as rast
FROM (SELECT rast FROM rasters.srtm1
    UNION ALL
    SELECT rast FROM rasters.srtm2) foo

And then used this function to create a slope for New York (this shapefile: https://data.cityofnewyork.us/City-Government/Borough-Boundaries/tqmj-j8zm)

WITH r AS ( -- union of filtered tiles
     SELECT
             ST_Transform(ST_Union(srtm.rast), 26918) AS rast
     FROM rasters.srtm as srtm
     JOIN borough_boundaries ny
             ON ST_DWithin(ST_Transform(srtm.rast::geometry, 26918), ST_Transform(ny.geom, 26918), 1000)
   ), cx AS ( -- custom extent
           SELECT
                   ST_AsRaster(ST_Transform(ny.geom, 26918), r.rast) AS rast
           FROM borough_boundaries ny
           CROSS JOIN r
   )
   SELECT
           ST_Clip(ST_Slope(r.rast, 1, cx.rast), ST_Transform(ny.geom, 26918)) AS rast
   FROM r
   CROSS JOIN cx
   CROSS JOIN borough_boundaries ny;

The query runs for a few minutes and produces some raster tiles, but when I try to visualize in QGIS I get the error that the raster contains different SRIDs and therefore it cannot load it. Is there a mistake in the query or a way to output raster tiles with the same SRID?

I'm using PostgreSQL 11 with PostGIS 2.5


texts: How to convert from per-hour to per-month units, respecting the variation in length of individual months?

texts: The GPM: Monthly Global Precipitation Measurement (GPM) v6 dataset consists of monthly images that have precipitation band with data in units of mm/hr. How can the precipitation data be converted into units of mm/month, respecting the variation in the length of individual months?


texts: Creating Raster that sums up input polygons using QGIS?

texts:
I am using QGIS and I have a shapefile with buildings. From these buildings, I have built Buffer with variegating radius. The Buffer also have data about the demand for Parking spots for each building.

How do I build a raster that accumulates the data of the Buffer?

Each pixel should sum up the attribute (demand of Parking spots) of all the buffer zones that overlay each other on this particular point(pixel).

My Aim with this raster is to show: "From this pixel I can serve so and so much parking demand (from the surrounding buildings)".


texts: PyQGIS (3.6.3 Noosa) Why is self.iface.mapCanvas().setExtent() being overridden in plugin?

texts:
I am writing a plugin that first adds the Google Imagery XYZ tile layer and then adds various other layers - some from a GeoPackage and others by creating their geometry in PyQGIS.

The last function of the plugin is to zoom to one of my created layers. However, no matter what I try, the map canvas shows the entire world (extent of the tile layer).

Printing the mapCanvas().extent() and the desired layer extent at the end of the plugin execute function return the desired extent but the map shows the tile layer extent.

I have tried without the tile layer and the canvas is set to the extent of the layer with the largest extent.

I have a push button in my plugin dockwidget that zooms to the layer extent using exactly the same functionality that is at the end of my plugin execute function. The button works fine after the plugin has finished processing.

# add the tile layer before the plugin does its processing
tile_url = 'type=xyz&url=https://mt1.google.com/vt/lyrs%3Ds%26x%3D%7Bx%7D%26y%3D%7By%7D%26z%3D%7Bz%7D&zmax=18&zmin=0&crs=EPSG3857'
rlayer = QgsRasterLayer(tile_url, 'Google imagery', 'wms')

QgsProject.instance().addMapLayer(rlayer)

# after a few seconds processing the canvas is supposed to zoom to the layer extent
# the print statements give the same extent but the canvas is showing the entire tile layer (or layer with largest extent)
# lyr is a vector layer being generated when the main function of the plugin is run.
# The following lines are at the end of the same function.

print(lyr.extent())
self.iface.mapCanvas().setExtent(lyr.extent())
self.iface.mapCanvas().refresh()
print(self.iface.mapCanvas().extent())

Is there a reason why the self.iface.mapCanvas().setExtent() is being overridden?


texts: Managing multiple Python interpreters (ArcGIS, OSGeo, etc.)

texts:
I currently have several Python interpreters installed through ArcGIS, OSGeo and Anaconda.

C:/OSGeo4W64/bin/python.exe
C:/Python27/ArcGIS10.4/python.exe
C:/Users/[UserName]/Anaconda3/python.exe

According to my understanding, each of these correspond to different versions of Python, and also have different libraries.

Consequently, if I want to run a script that uses arcpy, I use the Python interpreter. If I want to access the osgeo library, I use the OSGeo interpreter, and so forth.

This seems inefficient.

For instance, I only have one version of R which I update regularly and use for all my R libraries.

Could I link the Anaconda interpreter to the arcpy and osgeo libraries somehow or is there a reason why they require their own interpreters?


texts: Convert multiple JSON files to a single shapefile or geodatabase feature in ArcGIS or Python

texts:
I have been querying one feature in Esri ArcGIS REST Service, but because it is a large data, I had to query it several times (1000 features query limit). This results in many JSON files as parts of one that is supposed to be a single feature. Convert from JSON to Feature only allows one JSON. so

  Condition: I have many Esri JSON files

  How do I bulk-convert multiple JSON files and merge it into a single shapefile/geodatabase feature?

Is there a single way to do this using ModelBuilder in ArcGIS Desktop?

I am familiar with Python (GeoPandas) but not familiar with ArcPy.


texts: Adding new text field using if/then python statement?

texts:
I was hoping someone could help me with an issue I've been having while trying to aggregate data into a new column in an attribute table in arcgis. This is part of a model built using model builder. I added a new column, but the model always fails when it reaches the part where it has to populate the new column.

I'm trying to group specific land uses into a new land type (for example: field, pasture, and crop should all be coded as 'Agriculture' in the new column). I have been using the field calculator tool to populate the new column, but it's not working. I think either field calculator is the wrong tool to use in this scenario, or my python is wrong. Any advice? My python is below.

def landuse_groupFromName(LUDCODE):
if LUDCODE in ['23', '1', '36', '35', '2']:
return 'Agriculture'
elif LUDCODE in ['15', '16', '5']:
return 'Commercial/Industry'
elif LUDCODE in ['11', '10']:
return 'High_Density_Res'
elif LUDCODE in ['13', '38']:
return 'Low_Density_Res'
elif LUDCODE in ['12']:
return 'Med_Density_Res'
elif LUDCODE in ['34', '26', '29', '7', '8', '17', '9']:
return 'Recreation'
elif LUDCODE in ['18']:
return 'Transportation'
elif LUDCODE in ['31']:
return 'Urban'
elif LUDCODE in ['39', '19']:
return 'Waste'
elif LUDCODE in ['20']:
return 'Water'
elif LUDCODE in ['37', '14']:
return 'Wetlands'
elif LUDCODE in ['40', '3', '4', '6', '24', '25']:
return 'Wild'
else:
return 'undefined'


texts: Exporting subset of data but maintaining relationship class properties within new file geodatabase using ArcPy?

texts:
I have created a script tool that allows the user to choose a spatial extent and export several feature classes within that extent and (what is supposed to be) a related table to a new file geodatabase.  Within the master database where the exports are coming from, there is a relationship class between one of the feature classes and the table.  I am unsure of the best way to maintain the relationship (or create a new one) so that the records in the exported table reflect the related records that were deleted during the spatial extent clipping.

I have attempted to export the feature classes and table and then create a new relationship class within the geodatabase.  When I do this, the relationship works for any records deleted manually afterwards.  However, there are still a bunch of records in the table (destination table) that do not exist in the feature class (origin table) because they were not updated when the relationship class was created.  I have also tried creating the relationship class in the GDB prior to copying the features and table into it - this gives the same result.

Is there a specific way to make sure a relationship class works retroactively?  Or is there a better way to do this?

Here is my code for this part of the script:

output = arcpy.CreateFileGDB_management(workspace, gdb_name)
gdb = output.getOutput(0)

# make NHA_Core feature layer from selection and copy to new GDB
core_lyr = arcpy.MakeFeatureLayer_management(NHA_Core, "core_lyr")
core_selection = arcpy.SelectLayerByAttribute_management(core_lyr, "NEW_SELECTION", '"STATUS" = \'C\' or "STATUS" = \'NP\'')
if len(clip_lyr) > 0:
    core_selection = arcpy.SelectLayerByLocation_management(core_lyr, "INTERSECT", clip_lyr, "", "SUBSET_SELECTION", "")
else:
    core_selection = arcpy.SelectLayerByLocation_management(core_lyr, "INTERSECT", interactive_selection, "", "SUBSET_SELECTION", "")
Core = arcpy.CopyFeatures_management(core_selection, os.path.join(gdb, "NHA_Core"))
arcpy.DeleteField_management(Core, deleteFields)

Species = arcpy.CopyRows_management(NHA_SpeciesTable, os.path.join(gdb, "NHA_SpeciesTable"))

arcpy.CreateRelationshipClass_management(Core, Species, os.path.join(gdb, "NHA_Core_Species"),
"COMPOSITE", "NHA_SpeciesTable", "NHA_Core", "FORWARD", "ONE_TO_MANY",
"NONE", "NHA_JOIN_ID", "NHA_JOIN_ID")


texts: Add xy Data from csv adds data but changes xy values?

texts:
I have an issue with adding xy data from csv. files.
Data in csv properly formatted (name, x-point,y-point,z-point).

After importing them, when I click on Identify it shows location x and y coordinates that differ from the ones I tried to import.

Hopefully this picture will show what i am trying to say:

Any leads?

Ok. I open arcmap, set coordinate system, click on "Add xy data", point to csv file, select fields that contain x y and z coordinates, load it. Then I go to export data to shp file. Once the ship file is loaded, I click on one of the points to display xyz values. X and y are different then the values in csv file, z is ok


texts: Correlation between a nominal and a continous rasters in R

texts: I have two rasters of which one has nominal scale and the other has continuous scale. I used rasterCorrelation() function from spatialEco package in R and got a sensible map as well. But as one raster is at nominal scale, I am not sure about how do we explain it? rasterCorrelation() uses a matrix cum window of some specified scale. Statistically it's not possible to calculate correlation of nominal vs any variable because of zero variance. Then what is rasterCorrelation() function doing?


texts: Speeding up Raster to Polygon Conversion in ArcMap/arcpy

texts:
I am trying to convert a series of big binary rasters describing data in grey (value 0) and nodata in green (value 1) into polygons in order to derive the extent the data pixels are covering (see also this question: Creating index layer around rasters to rename files using ArcPy with ArcMap?)

This conversion process takes around half an hour for each raster dataset, probably because of the rather large file size and number of columns and rows. I enabled simplify polygons in the Raster to Polygon tool. The aim is to get out a raster data extent such as the blue sample shape in the screenshot below.

Is there a way to speed up this process or an alternative approach, e.g. by using numpy arrays? Or maybe a way to generalize the raster as I am sure the complexity of the raster leads to a longer raster to polygon conversion processing time.


texts: Understanding Google Maps overlay type from object (polygon or polyLine)?

texts:
var Line = new google.maps.Polyline({
    path:predectionArray,
    editable: false,
    strokeColor: '#FF0000',
    strokeOpacity: 1.0,
    strokeWeight: 2
  });

var Rect=new google.maps.Polygon({
    path:predectionArray2,
    editable: false,
    strokeColor: '#FF0000',
    strokeOpacity: 1.0,
    strokeWeight: 2
  });

the object Line is poly Line and Rect is Polygon .The question is how can i retrive it's type from Object name ,Like Rect.type or Line.type ?


texts: How to reorder output string in terms of values from a QTableWidget?

texts:
I have a QTableWidget which contains the names of layers and values. The code below creates a range of values from a list, iterates through some groups whilst creating a dictionary containing the layer name and value from the QTableWidget.

The output is a simple SQL expression which contains the layer name and its value in order of when the layer appears in the table. I want to change this order so that the layer with the higher values are included first.

Here is the code I used:

def func(self):
    root = QgsProject.instance().layerTreeRoot()
    group_name = root.findGroup('Group')
    qTable = self.dockwidget.tableWidget
    delim = ' WHEN '
    # Number of layers in 5 groups
    layer_count = [15, 14, 1, 2, 1]

    top = list(np.cumsum(layer_count))
    bottom = list(np.subtract(np.cumsum(layer_count),layer_count))
    pairs = zip(bottom,top)
    ranges = iter([range(x[0], x[1]) for x in pairs])

    # Iterate through the 5 groups
    for group in group_name.children():
        value_dict = {}
        for row in next(ranges):
            first_item = qTable.item(row, 0)
            first_name = str(first_item.text())
            second_item = qTable.item(row,1)
            value = str(second_item.text())
            value_dict[first_name] = value

        layers = [g.layerName() for g in group.children()]
        fscorestrs = {k:delim.join('"{layer}_Score" = {N} THEN "{layer}_Score"'.format(
                layer = l, N = value_dict[first_name]) for l in layers) for k in range(1)}
        formula = ("CASE WHEN " + fscorestrs[0] + " ELSE NULL END")

For example:

Here is a simple layout for the QTableWidget:

If I use the above code, the formula reads:

CASE
WHEN "Layer_1_Score" = 0.1 THEN "Layer_1_Score"
WHEN "Layer_2_Score" = 0.7 THEN "Layer_2_Score"
WHEN "Layer_3_Score" = 0.5 THEN "Layer_3_Score"
WHEN "Layer_4_Score" = 0.8 THEN "Layer_4_Score"
WHEN "Layer_5_Score" = 0.6 THEN "Layer_5_Score"
ELSE NULL
END

What I am looking for is to instead produce an output where the layers are ordered in terms of the value, so the above should look like:

CASE
WHEN "Layer_4_Score" = 0.8 THEN "Layer_4_Score"
WHEN "Layer_2_Score" = 0.7 THEN "Layer_2_Score"
WHEN "Layer_5_Score" = 0.6 THEN "Layer_5_Score"
WHEN "Layer_3_Score" = 0.5 THEN "Layer_3_Score"
WHEN "Layer_1_Score" = 0.1 THEN "Layer_1_Score"
ELSE NULL
END

How could I achieve this?


texts: Issues when downloading a WMS layer from Geoserver as KML

texts:
I am trying to download a Geoserver layer as KML using the KML Reflector. I actually use Leaflet and a plugin (leaflet-area-draw) in order to let the user design a rectangle and then pass the BB to the URL.
The URL request I make is WMS and it looks like this:

http://localhost:8080/geoserver/wms/kml?layers=my_workspace:my_layer&bbox=-124.73,24.96,-66.97,49.37

I face two issues. First the whole layer is downloaded although I restrict it with the BBox. I wonder if this is related cause its a WMS request and maybe I should try to make a WFS one. Is it possible to make a WFS request to get KML?

The second issue appears when the projetion of the layer in geoserver is other than EPSG:4326. Then I get nothing back. What is the right approach to this? Should I use a Leaflet plugin to convert the BBOX coordinates from EPSG:4326 to the layer's srs?

EDIT

I tried using KML scoring in order to get the data as vectors with the hope, that they will be rendered correctly and only the ones in BB will be shown. But it didn't work. I get all the features back.

http://localhost:8080/geoserver/wms/kml?layers=geonode:afg_wfp_acces_byair__byroad_20161801_1&mode=refresh&kmscore=100&bbox=68.0712890625,32.24997445586331,72.4658203125,35.88905007936091


texts: Condition and focal statisitics on multiple rasters

texts:
I'm trying to write a code in Python for several rasters, however I'm wondering if there is a faster way to do this.

I'm generating several rasters from a conditional statement for each value in the raster file. Then from these rasters I'm doing focal statistics.
I'm saving all the interim rasters from the conditional statement and using them in focal statistics. Is there a way to avoid that? Can I make the code below more automatic?

# Import arcpy module
import arcpy
from arcpy.sa import *
#
arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "D:\\_Work\\Test-Model.gdb"
#
CropMap = "RasterFile"
#
#
# Step 1
# Select crops
select_crop = "select_crop_136"
arcpy.gp.Con_sa(CropMap, "1", select_crop, "0", "Value = 136")
#
select_crop = "select_crop_146"
arcpy.gp.Con_sa(CropMap, "1", select_crop, "0", "Value = 146")
#
select_crop = "select_crop_147"
arcpy.gp.Con_sa(CropMap, "1", select_crop, "0", "Value = 147")
#
#
# Total Crop
tot_crop = "tot_crop_136"
arcpy.gp.Con_sa(CropMap, "1", tot_crop, "1", "Value = 136")
#
tot_crop = "tot_crop_146"
arcpy.gp.Con_sa(CropMap, "1", tot_crop, "1", "Value = 146")
#
tot_crop = "tot_crop_147"
arcpy.gp.Con_sa(CropMap, "1", tot_crop, "1", "Value = 147")
#
#
# Step 2
# Focal Statistics
#
neighborhood = NbrRectangle(1000, 1000, "MAP")
#
inRaster_sel_136 = "select_crop_136"
foc_select_136 = FocalStatistics(inRaster_sel_136, neighborhood, "SUM","")
#
inRaster_tot_136 = "tot_crop_136"
foc_tot_136 = FocalStatistics(inRaster_tot_136, neighborhood, "SUM", "")
#
# and so forth for all the rasters in Step 1


texts: Preventing loss of detail using "Clip Raster by Mask Layer" QGIS

texts:
I am using a raster file of satellite imagry as input for a 'Clip Raster by Mask Layer' process and a shape file polygon as my mask layer. Each time I run the process the output clipped layer is slightly brighter and slightly less detailed than the input.
I have checked the box for keep resolution of input raster so I was wondering if anyone has any ideas why this might be happening?
I'm using QGIS 3.12


texts: Computing area of raster cells with value "1" in QGIS?

texts:
I have two layers, a raster and a vector. The raster has values 0 and 1. The vector layer has an "Area" field for the total area of the shapefile/boundary. I want to compute the area covered of the raster valued with 1 against the vector layer.

Here's an image. What are the steps to achieve this?
The 1's are in white.

I am using QGIS 3.10 on Windows 10.


texts: Quick methods to make a shapefile/feature not m-aware?

texts:
Just wondering if anyone knows of a quick method to disable m-aware shapefiles/features? Searching Geonet, ArcDesktop help documents have proved unfruitful.

I'm currently doing copy features and changing it in environments. I could save a ton of time if I disable this in one go.

Anyone else have experience with this?


texts: Automatically align the border of one shapefile to overlap the border of another

texts: Can I automatically align the border of one shapefile to overlap the border of another, without editing each vertice individually?


texts: Plugin not working after building through Plugin Builder in QGIS

texts:
I have been trying to build a plugin through Plugin Builder 3 and been facing this issue then.

I have been facing the pyrcc5 issue but I was able to solve it through this link:
Fixing “The resource compiler pyrcc5 was not found”

This is the Python error:

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

ModuleNotFoundError: No module named 'split_layer.resources'
Traceback (most recent call last):
  File "C:/PROGRA-1/QGIS3~1.3/apps/qgis/./python\qgis\utils.py", line 334, in _startPlugin
    plugins [packageName] = package.classFactory (iface)
  File "C:/Users/pabla/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\split_layer\ _init_ .py", line 35, in classFactory
    from .split_layer import Split_Layer
  File "C:/PROGRA~1/QGIS3~1.8/apps/qgis/./pychon\qgis\utils.py", line 737, in _import
    mod = _builtin_import (name, globals, locals, fromlist, level)
  File "C:/Users/pabla/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\splic_layer\splic_layer.py", line 26, in
    from .resources import *
  File "C:/PROGRA~1/QGIS3~1.8/apps/qgis/./pychon\qgis\utils.py", line 737, in _import
    mod = _builtin_import (name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'split_layer.resources'

Can anyone help me how to solve this error?


texts: Quartiling a vector buffer in ArcMap 10.7

texts:
For reference, I found another Q&A with an answer to this same question for QGIS but I am looking for an ArcMap solution: Quartiling a vector buffer in QGIS 3.6.

I am trying to figure out if there is a way to cut the buffers into the needed cardinal quadrants. The goal is to then dissolve the four buffers for a point into one buffer that shows how the radius changes for each quadrant.

For context, I am working with NOAA hurricane data that has wind radii data values that extend in four directions: NE, SE, SW, and NE. Each value represents how far in each direction certain levels of wind speed extended at different points (every six hours). We create buffers out of those wind radii data to figure out where certain levels of wind speed intersect with different counties of the United States.

Our former ArcGIS analyst created circle buffers by taking the minimum of the four directional values: That is, if at a certain point a hurricane had a certain category of wind speed’s extending 25 miles NE, 50 miles, SE, 75 miles NE, and 70 miles NW she would create a simple circle buffer using the minimum value: a buffer with a radii of 25 miles (ArcPy code for this method and illustration at the end of this post).

What I would like to change is illustrated at the end of this post: I would like to create for partial circles with different radii that extend according to the four distances in the NOAA data. What ArcMap tools could I use to create the partial circled with different radii for a single point? I have done some research on my own and the only resource I have noted are multiple ring buffers, which, to my knowledge, cannot be made into partial circles.

Even a method to separate a circle buffer into four quadrants would be helpful. I think from there I could work on manipulating the radii of the separated vectors.

Sample ArcPy code for current approach:

arcpy.AddField_management("radii34","Buff_34","DOUBLE")

arcpy.CalculateField_management("pts","Buff_34","min_val( !NE50WR!, !SE50WR!, !SW50WR!, !NW50WR!)*1000", expression_type="PYTHON_9.3",

code_block='def min_val(a,b,c,d):\n if a == "-999":\n return 0\n else:\n return(min(a,b,c,d))')

arcpy.Buffer_analysis("pts","Hurricane_Buff34","Buff_34",line_side="FULL", line_end_type="ROUND",dissolve_option="NONE", dissolve_field="",method="PLANAR")


texts:
I'm looking for a routing engine written completely in native JavaScript or Nodejs, which is distributed with an open source license
So far the solutions I've found are implemented only in Java, C++ or Python.


texts: In ArcGIS 10.x, how do I select features of a point shapefile whose z values are above or below the elevations in a DEM?

texts: I have a multi-point shapefile representing earthquake hypocenters (containing depth values) and I have a DEM raster representing an underground surface (the top of the PreCambrian basement). I want to select only the points whose depth values are greater than the elevations in the DEM. Is there a way to do this?


texts: Headless QGIS: silent fail on broken import?

texts:
I am using a headless QGIS Processing script and found out that import errors (e.g. due to a missing library) lead to a silent fail. Neither does the object returned by processing.runalg contain info about an error nor is anything printed to the console.

This is not desirable behaviour and cost me quite some time to even realize that something went wrong.

To reproduce this problem put the two example scripts in a folder named scripts and run run_script.py in its parent folder.

scripts/my_script.py

print "Hello world. This is the processing script speaking."

scripts/my_broken_script.py

import unknown_library

print "Hello world. This is the broken processing script speaking."
print "Unfortunately you won't see this because the import above will fail silently."

run_script.py

import sys

# Prepare the environment
from qgis.core import *
from PyQt4.QtGui import QApplication

# ironic: we enable the GUI for headless operation, otherwise we get many errors like this:
# QPixmap: Cannot create a QPixmap when no GUI is being used
guiEnabled = True
qgis = QgsApplication([], guiEnabled)

# supply path to where is your qgis installed
qgis.setPrefixPath("/usr", True)

# load providers
qgis.initQgis()

# Prepare processing framework
sys.path.append('/usr/share/qgis/python/plugins')
import processing
from processing.core.Processing import Processing
from processing.script.ScriptUtils import *
Processing.initialize()

scripts = ScriptUtils.loadFromFolder('scripts')
name2script = {script.name: script for script in scripts}
print 'available custom scripts: %s' % (name2script.keys())

processing.runalg(name2script["my script"])
processing.runalg(name2script["my broken script"])

I noticed that if I run the script in QGIS (with GUI) a window pops up and tells me "no module named unknown_library" as expected. Is there any way to find out if processing.runalg was successful programmatically for headless scripts?


texts: Drawing polygon constant km around point in lat/lon coordinates in QGIS3.12

texts:
I have a set of points in lat/lon coordinates over Europe representing pixels from a satellite dataset. I know each pixel has size 10x10km.
I want to draw a square polygon around each single pixel with the size 10*10km. I have tried reprojecting and then drawing a square around it using the "Rectangles, ovals, diamonds" tool. But every time I get this way too big square drawn around the results (see the image below). I am not exactly sure how to fix or even tackle this issue.

I have tried the methods suggested in this question as well.


texts: Bringing WFS from Geoserver on top of Openlayers

texts:
I am trying to put a WFS layer on top of a openlayer's template I found in openlayers examples. This is the link for the untouched version. And this is the pen I am working on.
I am actually only changing the url defined in function extent and the center of coordinates to be zoomed. I think I do something wrong with getfeature link formation. Or maybe something about the coordinate systems.

<!DOCTYPE html>


texts: How properly calculate BBox for mocked GetMap request?

texts:
I am trying to imitate GetMap requests for GeoServer and have an application which generates these requests (requests are for cached layers). The problem is that sometimes I get geowebcache-cache-result:HIT and sometimes geowebcache-cache-result:MISS with geowebcache-miss-reason:request does not align to grid(s) 'XXXX'. I understand that for these "miss" requests I have generated invalid BBox parameter. For example (pay attention to marked numbers):

  ...REQUEST=GetMap&BBOX=55|6|865.8627022989,6104809.344416159,55|8|155.8627022989,6106099.344416159...

results in HIT

when

  ...REQUEST=GetMap&BBOX=55|5|865.8627022989,6104809.344416159,55|7|155.8627022989,6106099.344416159...

results in MISS with geowebcache-miss-reason:request does not align to grid(s) 'XXXX'.

What would be the formula for generating correct BBox?

NOTE: I have checked and the generated BBox coordinates fit inside layer's and grid extent.

UPDATE: Question also could be formulated as - how to validate if BBox coordinates are aligned to the grid that is used?


texts: Filtering layer by Intersects polygon on WFS 2.0.0 in Geoserver


texts:
Here is picture form QGIS.
There are three layers. Blue dotted lines are warehouse rack layer.
Green one is rack sections layer (routable network, created with pgr_createTopology).
And third layer is points (from _vertices_pgr table) with sections names.
It is saved as project. And the picture is generated with the command:

qgis-bin-g7.exe --snapshot print_points_wh.png --width 300 --height 1000 --project print_points.qgz --nologo --noversioncheck --noplugins

Layer three data source is table. When I need new plan, I just fill the table with required point ID`s. And then run command above.
The problem is that a GUI is required.

Then I tried to generate image with PostGIS tools.

select ST_AsPNG(ST_AsRaster((select ST_Union(the_geom) from print_points_geo), 400,1000))

Result is  Also these image dots are very small and only 3 points are visible. Two are missing.

If I add table print_points_geo to QGIS, result is

Is it possible to generate image like in first example from PostgreSQL without GUI tools?

Better would be to do it with PostGIS tools.
An option is to generate racks in one image, sections in another, and points in third image and then combine all images together. But I cannot figure out how to enlarge points in ST_AsRaster function.
Basically rack and section images are static.


texts: Is it possible to do directional variograms and kriging in SAGA or QGIS?

texts:
I'm trying to create directional variograms and do kriging for anisotropic soil data in SAGA or QGIS.  Aside from setting the direction to quadrant in SAGA, I can't seem to find a way to do this.  Is it possible to include a specific direction, set tolerance, etc.
I've looked at the manual and some tutorials but I can't find good examples of specifying the variogram direction.  I see that Surfer does have this functionality but I would like to use open-source software.


texts: QGIS Lost precision when converting raster to vector

texts: I am after running ordinary kriging, converting raster to vector I lose decimal value and only get whole number.  I have added first the raster then the vector image.


texts: Calculate values in one field based on values in another field

texts: I have one field with Zoning Codes (i.e.: RS) and I want to update another field with the string "Residential Single Family". The new field is currently NULL, I would like to add the full text string of the zoning codes. The field name containing the current code value is: [ZONING] and the new field is [ZONING_CLASS].


texts: Rasterizing lines significantly slower than polygons

texts:
I'm using GSHHG data to create land masks. Rasterizing the original polygons with GDAL_rasterize (via Python) is very quick. However, when I convert the polygons to lines and rasterize that, the process is much slower: seconds for the polygons vs minutes for lines.

Is there a particular reason for this e.g. some optimization I'm missing?

(in case you're wondering why I'm rasterizing the lines, it's to get a raster showing only cells with coastline in them)


texts: Graphical modeler in QGIS. I can't find a model

texts:
I'm a doing a tutorial for QGIS and at one point I have to use 'Join attributes by location' model but I do not find it in the Graphical modeler.

Does it have to be there?

I have to download it from somewhere?

I am using QGIS 2.8.1.


texts: Can I create a TIN from contours lines (shp) in QGIS?

texts:
I need to do a slope analysis, and would like to make a TIN from 2m contours lines.  Or do a slope analysis directly from contours? The only DEM available is not precise enough.  I do not have any elevation raster files.

And how do I then make slope classes?  The TIN properties does not seem to be able to classify anything other than elevation.


texts:
I am using OpenLayers to display a map and I have vectors layer with some features that I'm loading with AJAX from MySQL, and I'm using GeoServer to load all the background from shapefiles (like buildings, streets etc.) with default style, like this:

var generalLayers = new OpenLayers.Layer.WMS(
                "general", "http://XXX/geoserver/city/wms",
                { "transparent": "true",
                    LAYERS: 'city:Fence, city:Buildings, city:streets',
                    format: format,
                    tiled: true,
                    tilesOrigin: map.maxExtent.left + ',' + map.maxExtent.bottom
                },
                {
                    buffer: 0,
                    isBaseLayer: false,
                    displayOutsideMaxExtent: true,
                    yx: { 'EPSG:2039': false },
                    visibility: true
                }
            );

The features from MySQL are loading after 1.7 seconds, but the WMS layers are completely loading after 19.6 seconds!! (12 seconds after the features), all this time I have only the features without the background.

I have VPS so it's now because the server.  why it taking so much time? How can I improve it? Is it because I'm loading from shapefiles? if I will import everything (the shapefiles) to the db and I will load the background with MySQL plugin for GeoServer as WMS layer its will be faster?


texts: CARTO not styling 36 total representatives

texts:
I'm extremely new to CARTO and I'm still learning the ropes. I currently have this territory map where the colors represent a representative in a county.

However, it isn't shading the 36 total representatives. It only shades 10 then groups the other 26 in Others.

This is my current layer style setup in Builder:

All I want to be able to do is a have a color for every representative and not have the Others category. I read that I may have to do something with CartoCSS, but it didn't guide me in the right direction.

What's the best way to go about this? If there's a better way for me to visualize what I'm trying to accomplish, please tell me.

I'm running CARTO self-hosted.


texts: Finding distance between raster pixels and line features in R

texts:
I have a raster file and a SpatialLinesDataFrame in R. I need to find the distance between each pixel (at the center) and the nearest point in the line file. I believe this is the functional equivalent of the "near" function in ArcGIS. Is there a package or function that does this in R?

The raster is of a landscape with the year of deforestation as the cell value. The lines file is a series of roads- 34 total. I'm trying to determine the probability of deforestation around the roads using the distance of a point from the road as a variable.

The raster structure is:

# class       : RasterLayer
# dimensions  : 9003, 7906, 71177718  (nrow, ncol, ncell)
# resolution  : 0.00025, 0.00025  (x, y)
# extent      : -52.7825, -50.806, -1.125, 1.12575  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
# data source : U:\Jacy\Research\TS_line_analysis\GIS_Project\NE_Hansen
# names       : NE_Hansen
# values      : 0, 14  (min, max)


texts: I was wondering whether it is possible for a layered PDF created by MapInfo to hide the label of a layer in the layer list of the PDF, to prevent the user from altering the visibility of the layer, with the layer still being seen on the map?  I've heard you can right click the layer in the PDF document, and then select properties, but this seems to be disabled in my version of Adobe?


texts: Circle display on openlayers map (client-side vs server-side)

texts:
I'm developing a webGIS website where I draw circles from two sources:

a mapserver mapfile that "draws" them from PostGIS (red circles in image). The circles are created with ST_Buffer, eg

ST_Buffer(ST_GeomFromText('POINT(A B)', 4326), 0.072414686825054, 'quad_segs=32')

My mapfile serves in 4326 and 3857.

Via Openlayers (blue concentric circles in image), either via

var wgs84Sphere = new ol.Sphere(6378137);
new ol.Feature(ol.geom.Polygon.circular(wgs84Sphere, GPS_pos_4326, i * 5 * 1852, 64).transform('EPSG:4326', 'EPSG:3857'));

or via turf, then creating an OL polygon, by using turf's destination function.

I understand that projecting a circle in "Google Mercator" should produce ellipses when not on the equator, so the behavior is expected.

But how could I do the same thing via OL?


texts:
I need to make vector tiles but want the source to be like my data is.made up of polygons
Can GDAL ogr2ogr convert from polygon to line?


texts: Difference between OSM regional extracts unexpectedly large

texts:
I want to create a difference (change file) between two Europe extracts from OSM data. The files I want to compare are:

http://download.geofabrik.de/europe-180401.osm.pbf
http://download.geofabrik.de/europe-latest.osm.pbf

As the files are quite large, before running on them I have tried the processing on smaller extracts:

http://download.geofabrik.de/czech-republic-180401.osm.pbf
http://download.geofabrik.de/czech-republic-latest.osm.pbf

The workflow I use is:

first convert xxx-180401.osm.pbf to xxx-ref.o5m, then create the diff:

osmconvert europe-180401.osm.pbf -o=eu-ref.o5m

osmconvert eu-ref.o5m europe-latest.osm.pbf --diff -o=eu-changes.o5c

Or:

osmconvert czech-republic-180401.osm.pbf -o=cz-ref.o5m

osmconvert cz-ref.o5m czech-republic-latest.osm.pbf --diff -o=cz-changes.o5c

The process works fine with the Czech Republic data. The resulting o5c file looks reasonable, its size is ~30 MB (the Czech Republic input file is ~700 MB).

I get strange results with the Europe extract. The resuling o5c file is ~40 GB, while the Europe input is ~20 GB. When inspecting the file, I have found many instances of data which was not changed on OSM in this year at all, like way 539372444.

I have also tried doing the comparison using Osmium instead of Osmconvert, but the result was the same, the change file was huge.

Am I doing something wrong, or are the Europe extracts unsuitable for the comparison for some reason?


texts: Understanding "create-grid" output and getting grid points' coordinates

texts:
I created the following grid points by using the "create grid" pluggin from the processing toolbox, then "Intersection" to clip from a mask. My idea was to calculate the coordinates of these points with "Add X/Y fields to layer".
The issue is that I don't understand my attribute table. I was expecting to have 2 fields, but I have no idea what "top, right, left, bottom" represent. The issue is that because I have 4 fields, I fail to calculate the coordinates of the points with "Add X/Y fields to layer".
How could I calculate the coordinates of my points?
Why does the output of "Create grid" give 4 columns of data (aside from the indexes etc...)?


texts: Generalizing micro terrain in a slope raster

texts:
I am attempting create a map portraying avalanche risk levels: key to this is slope, aspect, elevation.

Slope output is too much of a checkerboard output, creating islands of low slope surrounded by high, etc. What is the best way to eliminate these "pockets"? Also, any tips on how to make it so any safe terrain that may fall below dangerous terrain is populate as dangerous?

I have some Python experience, but I am stumped as of now. I have used a 10m DEM thusfar.


texts: Selecting features using attributes of another file in QGIS

texts:
I have a table of attributes in which a field contains combinations of IDs from another shapefile (ex. '[6, 4]' where 6 and 4 are each individual IDs in my shapefile). Here's an example for visualization:

(the one in the left would be the aforementioned attribute table and the one in the right the shapefile).

I was wondering if there is some way I could select the features in my shapefile using the combinations in the separated table. Mostly I'd like to know if it's possible to Select by Expression using "what's selected in another table" as input.


texts: Setting default user role in ArcGIS Online to custom role?

texts:
I have created a custom role for ArcGIS online for all users in my organization.

Does anyone know how to set it as the default role when new users join via Enterprise login?

Right now the default role is "User".


texts: How can I count/find overlaps in ArcMap 10.4 but between area and line features (for example buildings and roads? Using topology I can find overlaps only between line-line or area-area feature.


texts:
I have a large point cloud in ASCII format that is aligned to the coordinate system I require. I also have a range of either .igs splines or reference surfaces .stp that I would like to filter the point cloud to in order to create new filtered point cloud groups with the name of the nearest reference object.

Is there a tool I can use in VTK, PCL, OpenGL etc that can parse the point cloud down rapidly? The cloud would be quite large at around 20gb and I would like to associate 20 reference features to the subsequent parsed cloud.

Would anyone have guidance for me to pursue this?


texts:
I have two layers, each contain hundreds of polygons. I'm trying to find a way to clip the data from one layer to each polygon of the other layer.
The problem is the clip features of the layer overlap. I would like the input layer to be clipped to each polygon separately and create duplicated in the overlapping areas. Is there an easy way of clipping the input layer to each polygon of the second layer separately? Unfortunately I'm not familiar with coding if that's an option.


texts: GPS coordinates are not matching up with street network shapefile despite having same projection

texts:
As some background, my colleagues and I used a shapefile of a city's streets as a reference to collect GPS data. They were written simply with a field for: street name, lat, and long. The many points make when looked at do seem to match to mirror the street network so I don't think the GPS is off.

However, I need to overlay the original street network shapefile with these new GPS points. The problem is that they are a bit far off from one another on the basemap despite having the same projected coordinate system. What I did was just adding the .csv file with the coordinates of each point to the map as XY data. According to these coordinates, the thousands of GPS data points are about a few hundred miles off from where the original street network layer is on the map. I'm sure this is an issue with projection, but I've been trying to troubleshoot it for an hour and haven't managed to get them to line up or even move.

I have checked the source and they are both the same. I have tried batch projecting them as well, but it just results in duplicates being placed upon each of the shapefiles.

Here are the contents of the Data Source of each layer.
GPS:

Data Type:  Shapefile Feature Class
Shapefile: path\to\data
Geometry Type:  Point
Coordinates have Z values:  No
Coordinates have measures:  No

Projected Coordinate System:    NAD_1983_StatePlane_Illinois_East_FIPS_1201_Feet
Projection: Transverse_Mercator
False_Easting:  984250.00000000
False_Northing: 0.00000000
Central_Meridian:   -88.33333333
Scale_Factor:   0.99997500
Latitude_Of_Origin: 36.66666667
Linear Unit:    Foot_US

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

Street Network

Data Type:  Shapefile Feature Class
Shapefile:  \path\to\file
Geometry Type:  Line
Coordinates have Z values:  No
Coordinates have measures:  No

Projected Coordinate System:    NAD_1983_StatePlane_Illinois_East_FIPS_1201_Feet
Projection: Transverse_Mercator
False_Easting:  984250.00000000
False_Northing: 0.00000000
Central_Meridian:   -88.33333333
Scale_Factor:   0.99997500
Latitude_Of_Origin: 36.66666667
Linear Unit:    Foot_US

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

I'm pretty new to GIS, so I'm hoping it may be something simple. Does anyone know of any way I could overlay these two layers? Is there something I'm overlooking?


texts: proj/pyproj converted coordinates offset

texts:
I'm using pyproj to convert some coordinates from WGS84 to OSGB36. When mapped out the converted and unconverted don't map to the same places on a map. The OSGB36 coordinates are mapped about 100 meters north west of where they ought to be.

Does anyone know what might be wrong here?

Here's an example of what I'm doing in code

import pyproy
wgs84=pyproj.Proj("+init=EPSG:4326")
osgb36=pyproj.Proj("+init=EPSG:27700")
row = pyproj.transform(wgs84, osgb36, lng, lat)


texts: How to Import a TrueType Font (or OpenType Font) into a custom ESRI Style without doing it one symbol at a time?

texts: Creating a personal TrueType or OpenType Font with font creation software is easy. Pulling that creation into a custom ESRI Style as (e.g.) Character Marker Symbols is a tedious affair. Is there a tool or scripting method of populating a custom style with my new TTF of approximately 140 characters?


texts: QGIS - formula works in the attribute table fine but getting error when using the model builder

texts:
So I am using QGIS 2.18.15

Strange thing:

If I go via the attribute table into the field calculator and use following formula it works just fine:

  if( "DB_arit" IS NULL, 0, "DB_arit" )

Explanation: It searches for NULL values in the column "DB_arit" and replaces them with "0", all other values stay the same.

Now if I use the modelbuilder I get the following error:

  Algorithmus Test Field calculator NULL to 0 in DB_arit startet...
  Bereite Algorithmus vor: QGISFIELDCALCULATOR_1 Verarbeite Field
  calculator Field calculator NULL to 0 in DB_arit [1/1] Parameters:
  INPUT_LAYER =replacedbyme.shp, FIELD_NAME =DB_arit, FIELD_TYPE =0, FIELD_LENGTH
  =10, FIELD_PRECISION =2, NEW_FIELD =False, FORMULA =if("DB_arit" IS NULL, 0, "DB_arit" ) Fehlgeschlagen

  Fehler in Algorithmusausführung Field calculator Field calculator NULL
  to 0 in DB_arit QgsFeature.setGeometry(QgsGeometry): argument 1 has
  unexpected type 'NoneType' See log for more details

(Don't wonder aboout "Input_Layer=", I edited it just for this post)

Can anyone explain or help?