texts: How to ignore "handle bad layers" in QGIS?

texts:
I have a QGIS (2.14.3-Essen) project that I use both at work and at home. I keep all but two layers with the project.

The two remaining layers are very large - a DEM and hillshade totaling more than 20 GB - and they are stored in different paths at work vs. home.

I would like to ignore the Handle Bad Layers dialogue and keep the broken paths for those layers. I often don't care they are bad; I want to make changes to other parts of the project.

Is there a way to do so, to ignore or skip the Handle Bad Layers dialogue so I don't lose the reference to the layers if I don't update them?


texts: Calculating distances from one start point to one end point using batch data in MapInfo?

texts:
Im trying to calculate a batch of data with a pupil unique identifier, their home postcode and their adoption placement postcode.

I would like to find the distance between the home and the placement of the pupil. I do have RouteFinder but only a small road network so I'm happy to just get a straight line distance measurement.

I have tried creating two sperate tables with the geocoded postcodes, one for home postcodes with the unique pupil ID and one with their placement postcodes with the ID.

Has the distance calculator got a capability of just calculating the home postcode point to the placement postcode point where the identifiers match.

I know i can get every possible combination and then filter it out in Excel using a formula however i can only do this for small batches of info.

Can anyone let me know if there is a better way to do this with a large data set?


texts:
I'm trying to create a one-way replica between two ArcGIS Enterprise Geodatabases. The process works when I select only the tables in the parent. The replica will successfully copy all datasets to the child database. However, when I select the feature classes (there are hundreds of them), the process fails after about 17-18 minutes or so.

Replica Settings:

Error Message:

  ERROR 000582: Error occurred during execution. Item with the same path
  name already exists.

The child database is essentially empty and doesn't contain any feature classes yet. My first thought was that the process was truncating long feature class names, however, the feature class name length limit is 255 characters and the longest one (including database name and schema is 95 characters).

Both databases are on different servers. Both servers are are Windows Server 2016. Databases are SQL Server 2017.

I had a look at this site, where one suggestion is to create a replica on each individual dataset to see which one is causing the error. However, there are over 600 feature classes. Is there an easy way to figure out what is going on here?


texts: Layer not visible when rendering from Python script in QGIS

texts:
I am running a script from the Python console in QGIS. I have stripped the script down to highlight the problem. Basically, when I render and save an image of a shapefile that is displayed with points they do not appear in the saved image. They will appear in the map canvas (see image below). Interestingly, if I have a vector that I fill (e.g., administrative boundaries) they will be visible in rendered image. Here is a screenshot of the map canvas after running script:

and here is my code:

from PyQt4.QtCore import *
from PyQt4.QtGui import *

#input Layer
lyrPts = QgsVectorLayer("SomeFile.shp", "Qflu", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(lyrPts)

# create image
img = QImage(QSize(8000, 6000), QImage.Format_ARGB32_Premultiplied)

# set image's background color
color = QColor(255, 255, 255)
img.fill(color.rgb())

# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)
render = QgsMapRenderer()

# set layer set
layers = QgsMapLayerRegistry.instance().mapLayers()
lst = layers.keys()
render.setLayerSet(lst)

# set extent
rect = QgsRectangle(lyrPts.extent())
rect.scale(1.1)
render.setExtent(rect)

# set output size
render.setOutputSize(img.size(), img.logicalDpiX())

# do the rendering
render.render(p)

#end rendering
p.end()

# save image
img.save("image.png","png")


texts:
Apparently, the new released SNAP 6.0 allows the sen2cor plugin to be installed and run directly within the SNAP. Indeed, SNAP installs the bundle when try to run it.

However, it doesn't work for me. Although a new folder is generated (S2A_MSIL1A...from the S2A_MSIL1C), the folder does not contain the corrected bands.

Is there anyone that succeeded with it without using Anaconda nor the standalone sen2cor tool?

I enter the image through the MTD_MSIL1C file within the download folder S2A_MSIL1C_...

Then I run the sen2cor  and a new folder is generated named S2A_MSIL2A_...

But inside the latter folder I just have two empty folders: AUX_DATA and DATASTRIP

Any suggestion?


texts: Cannot access MODIS Download from FTP using Python

texts:
I tried to access the FTP site of MODIS data using Python ftplib module but always fail to connect. I am connecting using: ftp = ftplib.FTP("e4ftl01.cr.usgs.gov")

here is the error: A connection attempt failed because the connected party did not properly respond after a period of time, or establishment connection failed because connected host has failed to respond.

However, if I access the ftp site thru my web browser, I was able to connect and download files.


texts: Why is Earth Engine Export.map.toCloudStorage taking much longer than before?

texts:
I've created a script to identify the nth least cloudy mosaic for Sentinel-2 data for a given place and time. A few months ago, I was able to execute an Export.map.toCloudStorage for a visualization of Sentinel-2 data for an area a little over half as big in less than 20 minutes with the same maximum zoom. However, this export (which is a rectangle extent about 40,000 sq km) has been running for 17 hours. It appears to be working,  as I can see some map tiles in the bucket, but just taking a really long time. Have I done something in this script that would cause that or has there been a change to Earth Engine services?

The script:

/* This script will export the nth least cloudy Sentinel-2 mosaics for a given
region of interest (ROI), year, and range of months* as map tiles to google
cloud storage */

// Select year

var year = 2019

// Select start month

var startmonth = 3

// Select end month

var endmonth = 6

// Load the ROI table
var table = ee.FeatureCollection("users/theronmb/TTCS_Area_Extent_for_Sentinel-2");

// Find the MGRS_TILES that overlap with the ROI for the time period of interest
var collection = ee.ImageCollection("COPERNICUS/S2")
  .filter(ee.Filter.calendarRange(year,year,'year'))
  .filter(ee.Filter.calendarRange(startmonth,endmonth,'month'))
  .filterBounds(table)
;

//print(collection)

// Get the unique MGRS tile list
var tiles = ee.List(collection.distinct('MGRS_TILE')
  .aggregate_array('MGRS_TILE'));

var orbits = ee.List(collection.distinct('SENSING_ORBIT_NUMBER')
  .aggregate_array('SENSING_ORBIT_NUMBER'));

/* Function to iterate through the image collection to pull out the nth
least cloudy granuales for each target tile.*/
var fun11 = function(orbit) {
  orbit = orbit;
  var fun1 = tiles.map(function(tile) {
    var S2 = collection
    .filterMetadata('SENSING_ORBIT_NUMBER', 'equals', orbit)
    .filterMetadata('MGRS_TILE', 'equals', tile)
    .sort('CLOUDY_PIXEL_PERCENTAGE', true);
    var mylist= S2.toList(3);
    var s2 = ee.Algorithms.If(ee.Number(S2.size()).gt(1),
      ee.Image(mylist.get(0)), //change to 1 or 2 for 2nd or 3rd least cloudy
      ee.Image(0));
    return s2;
  });
  return fun1;
};

// Apply function to orbits list
var list1 = orbits.map(fun11).flatten();

// Create image collection from image list and remove blanks
var collection1 = ee.ImageCollection(list1)
    .map(function(image) {
      return image.set('count', image.bandNames().length())
    })
    .filter(ee.Filter.eq('count', 16))
    .sort('CLOUDY_PIXEL_PERCENTAGE', false);
print(collection1)

var vizParams1 = {
  bands: ['B8', 'B3', 'B2'],
  min: [640, 524, 671],
  max: [4347, 1189, 1220],
  gamma: [1, 1, 1]
};

var vizParams2 = {
  bands: ['B8', 'B3', 'B2'],
  min: [1400, 276, 157],
  max: [3910, 1624, 1183],
  gamma: [1, 1, 1]
};

Map.centerObject(table, 8);
Map.addLayer(collection1, vizParams1, 'least cloudy', 0);

// Create visualization layers.
var sent1RGB = collection1.select(['B8', 'B3', 'B2'])
  .map(function(image){return image.clip(table)})
  .mosaic()
  .visualize(vizParams1);

Map.addLayer(sent1RGB);

// Export the image to Cloud Storage.
Export.map.toCloudStorage({
  image: sent1RGB,
  description: 's2_2019_lc1',
  bucket: 'ttcs',
  maxZoom: 15,
  minZoom: 8,
  skipEmptyTiles: true,
  region: table
});

*****Update*****

The export completed in 18 hours. I decided to test with a much smaller area (just 2600 sq km) and the script I have used in the past to export with quick results. It has already been running for 40 minutes. It appears that something has changed with regards to the resources that Earth Engine devotes to map exports. Also, even though I specified in the export task that the map tiles objects should be public in the export, they end up being private in the bucket, where as before they were correctly set to public.


texts:
Original question: I would like to sort my points in QGIS. See image

Restated question: The grid tool currently numbers by starting at R1C1 and going to R1Cx, then R2C1 to R2Cx. How do I number the points from R1C1 to R1Cx, then R2Cx to R2C1? Can the solution be applied after a regular grid has been destroyed (points removed) or does it need to be applied to a whole grid, assuming no option for this ordering scheme exists in the grid creation tool itself?

The point is to avoid the long 'row reset' or 'back to start of row' line as this wastes time and fuel, rather than moving directly down to the next row.


texts: Slight translation with coordinate-transformation from Gauss-Kruger to WGS 84. Possible Reasons?

texts:
I have a number of points in Gauss-Kruger, Zone 3 that I want to convert to display them within Google-Earth. Assuming WGS84 for Google-Earth the transformed coordinated differ from the original ones about 140m to the north-east.

I am doing the coordinate transformation using the "PROJ.4 - Cartographic Projections Library" (http://trac.osgeo.org/proj/) with the following proj4 strings:

Gauß Krüger, Zone 3: +proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs
WGS84:  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs - "

Could somebody check if my proj4 strings are correct (especially WGS84 for Google Earth) or suggest other causes of the error?


texts: display WMS layer from geoserver in mobile app with openlayers mobile

texts:
i try to display WMS layer from geoserver with openLayers Mobile,and its working great in the navigator (chrome)  but its not appeering in the simulator device (android)...
 & with another WMS layer like open street map , its working perfectly in both, web navigator and android device , so i can't figure out where the problem is , its been more then 10 days and i can't handle this problem..

this is my code:

var map;
    OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";

    function init() {
        map = new OpenLayers.Map('map', {
                controls: [
                    new OpenLayers.Control.Navigation(),
                    new OpenLayers.Control.PanZoomBar(),
                    new OpenLayers.Control.LayerSwitcher({'ascending':false}),
                    new OpenLayers.Control.ScaleLine(),
                    new OpenLayers.Control.MousePosition()
                ]
            });

            var wms = new OpenLayers.Layer.WMS(
                "OpenLayers WMS", "http://vmap0.tiles.osgeo.org/wms/vmap0",
                {layers: "basic"}
            );

           var osgeo = new OpenLayers.Layer.WMS("OpenLayers WMS",
                "http://vmap0.tiles.osgeo.org/wms/vmap0",
                {layers: 'basic'},
                {isBaseLayer: true, transitionEffect: 'resize'}
            );

          com = new OpenLayers.Layer.WMS(
            "communes",
            "http://localhost:8080/geoserver/wms",
            {layers: 'cite:matcommune_view', format: 'image/gif',transparent: true}
        );

       map.addLayers([osgeo, com]);

map.setCenter(new OpenLayers.LonLat(-104, 42), 3);


texts:
I have created a table join between a vector shapefile and a CSV file using pyQgs.

Now I want to edit the CSV file (not from inside QGIS), but when I double click on the CSV file, it says that the file is locked for editing by another process (which means by QGIS).

Is there any way I can edit the CSV file without losing QGIS or removing the CSV file from my QGIS project?

Once I edit the CSV, I just want to reload the edited file into QGIS project by another Python script.  I just do not want to close QGIS at any point.


texts:
I need custom rotation in feature text labels. Although this functionality is available in OL3, but at this moment I am unable to migrate to OL3 from Ol2.

I have tried changing SVG code at some levels, but nothing worked.


texts: Converting from ESRI Shapefile to Aircom Asset data format

texts:
I've shapefile contains vector layer of 3D buildings with their Z values.

I'm looking for help to convert this file to AIRCOM International’s ENTERPRISE format. The aim is to use converted vector data in Asset planning tool

Any tool or programmatic/Python solution is welcome


texts:
I am doing a undergraduate GIS project. It entails a paper route for newspapers in the Pittsburgh Area, so far I have plotted all of the points in a point layer. I have now plotted my points, what should be my next steps? I downloaded basemaps of the roads, and the data set for the pgh roads. I want to: find the shortest distance (elucian distance), I need to start a new network data set from what? (how)?

Please suggest anything else that would be useful or make this project better. Or a list or steps necessary for completing the project.


texts: Overlay a shape into a gridded shape and determine the percentage of each cell that are overlapped

texts:
I would like to know how I can overlay a shapefile (A) into another gridded shapefile (B) and extract the information about which shapefile (A) polygons falls in each cells (B).
And, more specifically, I would like to know how I can determine the percentage of each cell that are overlapped by a shapefile.


texts: Drawing polygon border on inside only in Mapnik?

texts:
I must draw a boundary of a natural reserve area in mapnik so that it looks like this:

It is a partially opaque green polygon with a wide line border (also partially opaque) on the inside of the said polygon.

Now, the problem is - how do I tell mapnik (in style XML) to draw the line only on the inside(!) of the polygon.

Mapnik wants to draw the line axially centered at the polygon boundary, so that half of it always "hangs" outside of the polygon like this, and it is no Ok:

The code for the result you see is this:

  <Style name="nature_reserve">
    <Rule>
        <PolygonSymbolizer fill="#4eaa2f" fill-opacity=".3" gamma=".2" />
        <LineSymbolizer stroke="#4eaa2f" stroke-opacity=".5" stroke-width="20" />
    </Rule>
  </Style>

  <Layer name="nature_reserve" status="on" srs="+proj=longlat +datum=WGS84 +no_defs">
    <StyleName>nature_reserve</StyleName>
    <Datasource>
      <Parameter name="type">shape</Parameter>
      <Parameter name="file">data/osm/nature_reserve</Parameter>
    </Datasource>
  </Layer>

I am a newbie to Mapnik, so that is about all I can think of. Any better suggestions?

Never mind the jagged lines, that is irrelevant now. All  I care about is getting the outline inside the polygon.

Can it be offset somehow?

I did not find any options for this in documentation.

They are doing it in OpenStreetMap, so there must be the way:


texts:
I am having problems trying to write ArcGIS Online Feature Service in update mode.

Our process is very simple. We published data from Oracle SDE to ArcGIS Online Feature Service with ArcMap. The source is in Place State projection, we set the data frame to WGS84-EPSG 4326 before publishing the service to AGO. We need FME to be able to update the existing AGO Feature Service when a change occurs in the corresponding SDE source layer (records number change). We want to keep the process simple for now, all I need is to write the existing Feature Service with the update and keep the WGS84.

I tried in two ways:

Read the SDE layer - Write existing Feature Service in ArcGIS Online (writer mode: insert)
Read both SDE feature and his corresponding Feature Service, use a FeatureMerger Transformer to check for updates and write the existing Feature Service in update mode (attached)

..but it doesn't work.


texts: Change No Data Value Geotif File with Qgis or Gdal

texts:
I have a little question for Change No Data Value Geotiff File with Qgis or Gdal by another value. Usually one imports or exports a file with no data value -9999, but sometimes we need to change that number for another on raster calculator the expression != but it seems not to work with No data value
any idea

Nice discussion, Carlos  but at same Kurt this option is for a range of data.
I was user syntax band math of ENVI and I thought to apply something like this on
Raster Calculator

( "MCD10A1_2018266@1" != -9999) * -8888 AND  ( "MCD10A1_2018266@1" >=  0 )  * "MCD10A1_2018266@1"

But I do not understand does not work!!


texts: Where can I find a tutorial for installing an addon to GRASS GIS that a six-year-old could understand?

texts:
I want to use the r.stream.distance addon in GRASS GIS 6.4 on ubuntu 10.04 (64 bit).

I've installed grass from the repositories and that works fine but I can't seem to figure out how to install this addon (or any other for that matter).

I've tried to run g.extension but I get a "command not found" error.  I've tried to install the addon from svn but the instructions are less than clear to me.

I've read through the documentation on the grass wiki and it seems to assume a lot of knowledge that I seem to lack.

So, does anyone know of a tutorial where I can walk through the steps (baby steps) of installing an addon in grass?

Thanks.


texts: stratigraphic ranking with field calculator?

texts:
I would like to populate a feature class float attribute with a numerical ranking of a text attribute of a geologic map. I have a field named [sheet_unit] with strings such as "H10, Q" and "J10, Pl" with I would like to relabel in another attribute [Stratigraphic_Rank] with a numerical value such as 0 and 10, etc. I was hoping to use an IF THEN statement but I cannot for the life of figure out how to use the field calculator to do this. Anyone willing to help me start?

This and my several variations of it don't work:

If [sheet_unit] = "H10, Q" THEN
[Stratigraphic_Rank] = 0
ElseIf [sheet_unit] = "J10, Pl"
[Stratigraphic_Rank] = 10
EndIf


texts: Weekly raster data for temperature and precipitation in US

texts: I am looking for weekly raster data for US for temperature and precipitation. All the major datasets PRISM, WorldClim etc. are either on a monthly scale or daily scale. I was wondering if anyone has aggregated data for weekly values from daily data or can point me to links where I can get it. I am looking for data for years 2005-2011. Thanks.


texts: Finding all the values within a band in google earth engine

texts:
I have an image composite which includes a day of year band. When I display the band using Map.addLayer and stretch by 100% in the Visualization Parameters I can see that the range of values is 151-281.

How can I get a list of all of the values within this range, i.e. a list of all the unique values of 'doy'? It should not be a long list as the composite image was derived from Sentinel 2 images restricted to a 5 month window. I just need to print the list. 'aoi' is a multi-polygon shapefile.

//Function to mask clouds using the Sentinel-2 QA band
function maskS2clouds(image) {
  var qa = image.select('QA60');

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

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

  return ee.Image(image.updateMask(mask)
  .divide(10000).copyProperties(qa)).set('system:time_start', qa.get('system:time_start'));
}

// Load Sentinel-2 TOA reflectance data.
  var S2 = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2017-06-01', '2017-10-15')
                  .map(function(image){return image.clip(aoi)})//To clip all images to aoi
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  //Select required bands only
                  .select('B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'QA60')
                  //Apply cloud mask
                  .map(maskS2clouds);

//Function to calculate redEdgeNDVI
var add_reNDVI = function(image){
  //Create band variables
  var redEdge = image.select('B5');
  var NIR = image.select('B8');
  var redEdgeNDVI = NIR.subtract(redEdge).divide(NIR.add(redEdge)).rename('reNDVI');
    return image.addBands(redEdgeNDVI);
};

//Map sentinel 2 collection to reNDVI
var S2_reNDVI = S2.map(add_reNDVI);

//Function to create DOY band
var addDate = function(image){
  var doy = image.date().getRelative('day', 'year');
  var doyBand = ee.Image.constant(doy).uint16().rename('doy');
  doyBand = doyBand.updateMask(image.select('reNDVI').mask());
    return image.addBands(doyBand);
};

//Map sentinel 2 collection to DOY function
var S2Date = S2_reNDVI.map(addDate);

// Make a "max redEdgeNDVI" pixel composite taking max reNDVI from all images
var S2_max_reNDVI = S2Date.qualityMosaic('reNDVI');

//Add DOY layer
Map.addLayer(S2_max_reNDVI, {bands:['doy'], palette: ['ffeda0', 'feb24c', 'f03b20']}, 'DOY for Max reNDVI pixel composite');```


texts: Where can I get the Cadinput plugin for QGIS

texts: From Stackexchange, I got to know that Cadinput plugin could do this work. I need to create a line feature specifying the desired length and bearing in an open source GIS (Preferably QGIS) as in CAD drawing. This is for the creation of floor plan of a building.


texts: Postgres casting between varchar and double precision suddenly stopped working

texts:
I have been importing data into PostgreSQL tables all day using ogr2ogr and then spatially enabling the data with no issues up to now. I have been cutting and pasting the code and this is what I have been using:

ALTER TABLE mytable ADD COLUMN geom geometry(POINT,4326);
UPDATE mytable SET geom = ST_SetSRID(ST_MakePoint(cast (lon AS double precision),cast (lat AS double precision ) ),4326);

On the last table, I suddenly get this error:

ERROR:  invalid input syntax for type double precision: "#VALUE!"

Nothing has changed, at least as far as I know, I was importing a csv file, same as all the others.

I checked the csv file for spaces or unusual characters but it looks in order

    Column lon, lon are of type character varying

here is a sample of the data

id,name,name_alt,lat,lon,description
1,Abbott Seamount,,31.8,174.3,A seamount located north of Midway Islands (US) in the North Pacific Ocean; name origin unknown.
2,Abra Canyon,,17.51667,120.36667,An undersea canyon located off the northwest coast of the Philippine island of Luzon in the South China Sea; named after the Abra River
3,Abraham Canyon,Abraham Sea Valley,52.61667,172.33333,An undersea canyon in the area of the Aleutian Islands in the North Pacific Ocean; named in association with Abraham Bay; Attu Island; Alaska (US).
4,Abrolhos Bank,,-18.5,38.75,A bank located off the coast of Brazil in the Atlantic Ocean; upon which sits the Arquipélago dos Abrolhos (archipelago).
5,Abrolhos Seamounts,Abrolhos Ridge; Cadeia dos Abrolhos,-17.33333,36.5,Seamounts located off the coast of Brazil in the Atlantic Ocean; named in association with the Abrolhos Bank and Arquipélago dos Abrolhos (archipelago).
6,Abubacer Ridge,,36.8,1.65,An undersea ridge located off the coast of Spain in the western Mediterranean Sea; named for Ibn Tufayl (known as Abubacer to Europeans); 12th-century Spanish-Arab philosopher and physician.
7,Aby Canyon,,4.25,-3.66667,An undersea canyon located off the coast of Côte d'Ivoire; named in association with Lagune Aby; Côte d'Ivoire (lagoon)

What am I missing here?


texts: Making virtual raster with average values

texts:
I'm working with many raster files. Is there a way to make virtual raster with average values in each pixel. (QGIS)

It's for an NDVI index of summer months over 5 years.


texts: how to make decent sense from scattered river depth measurements

texts:
I sampled river depth at scattered positions and I want to be able to make a map and GIS layers out of it.

My problems are:

As the river meanders the interpolation I do in scientific packages i,e. matlab generates values in points where there is no river therefore when I create a KML file out of it, it covers land areas. Does QGIS/GRASS have a better solution for that?
My GPS points accuracy of ~3m 90% of the time. Do I need to do any filtering for location?


texts: How to assign a unique ID to groups of items in the attribute table in ArcMap?

texts:
Using arcmap's field calculator I need to programatically loop through all of the entries in the field and for each name it gets its own unique id.  So names that are the same would get the same number. This is done through model builder.

def NameReclass(name):
    foreach name in species .... etc.

I'm just really confused how to do this without if/elif statements and specifying the string name == 1, etc. because it needs to be an id for each name shown.


texts: Unable to disable all layers visibility in GeoServer except admin

texts: I need to disable all layers visibility in GeoServer except admin. It works fine in GeoServer Webpanel but not reflect the same in ArcMap. Even the WMS links gives "Connection Error" when I change the ..r to admin


texts: Creating distinct polygon zones from interpolated raster in ArcGIS for Desktop?

texts:
I would like to create smooth/distinct polygon zones from raster data that was interpolated from sparse point data.

Is the only way to do this is to have a very small output cell size during the interpolation process?

I don't like the look of the raster data when I zoom in and see the grids from the raster data.

Is there a way I can eliminate these raster grids and have truly distinct polygon zones?

In addition, it would be nice if I can change the number of classes that make up the zones easily by simply changing the symbology and changing the number of classes used or even the classification method used.


texts: Calculating Stream Order Classification for NHD High Resolution vector data

texts: So I am trying to calculate stream order using the Strahler method. I have vector data from NHD in high resolution and was wondering if anyone has ideas on how to assign stream order using vector data. I prefer opensource solution. Any ideas?


texts: Obtaining road distances between each pair of points in point feature class using ArcGIS for Desktop?

texts:
In Calculating attribute for network distance between multiple points in ArcGIS Desktop?, an explanation has already been given how to find road distances from a certain point to each point in another feature class. I am now trying to find out how to find such road (not Euclidean) distance between each pair of points in the same (large) feature class. In other words, I am trying to obtain a matrix showing the road distance from each point in a feature class to all others.

What is the easiest way of doing this in ArcGIS?

I presume that this, too, involves Network Analyst, but I could be wrong about this.


texts: I would like to extract individual trees in ArcMap 10.1. I have 1.3 meter resolution DSM data created from aerial images. I used Crown Height Model to extract tree heights (as attached file), but I am stuck here. Is there any possibility to somehow extract highest points (single point for each tree) or delineate each tree based on their heights.


texts:
Desktop 10.3.1
GCS: WGS 1984

Created a Buffer of points based on a pre-specified radius from the points attribute table.

It appears all of the buffers were created with a decimal degree unit of measure.... How do I change the units of the buffer layer to miles?


texts: Assigning list or array to IRowBuffer Row without Row.set_Value()?

texts:
Is it possible to assign a list or array to IRowBuffer Row?

I have something like the following:
I have a fields (list) with values and i have to fill a Row with that values.

public void fillRow(IRowBuffer Row)
 {
  for (int i=0;i < fields.length;i++)
    {
      Row.set_Value(i, fields[i].toString()
    }
 }

Is there no way where i can just assign the fields to Row, without looping. What if my fields has like 500 elements and i have to loop 100x500? I know that Row has 2 requiredFields (ID, SHAPE)


texts: Save a shapefile with the original name as the imported file in R

texts:
I want to save as shapefile an object but I want the name of the shp to be the same (original) as the file I imported at the beginning of my code. I don't want to do it manually because at a later stage I will import hundreds of texts and I will want to save them as shp with their original names. Here is my code:

setwd("mydir")
d = read.table("vourlionas.txt",
               header = FALSE,
               sep="",
               col.names=c("x", "y"),
               fill=FALSE,
               strip.white=TRUE,
               stringsAsFactors = FALSE)

some data manipulation here

writeOGR(obj=test1, dsn="mydir", layer = "ORIGINALFILENAME"
               driver="ESRI Shapefile", overwrite_layer = TRUE)

I believe it has to do something with layer parameter. Any thoughts?


texts: How to tell SharpMap to display polygons based on z-order?

texts: I have a series of polygons which overlap in spatialite. Each has a an attribute providing it's z-order as well as the output color for that polygon. What do I need to do, processing-wise, to allow SharpMap to display the polygons based on z-order, along with the proper color for each polygon, all in a SINGLE layer?


texts: How to truncate the %Name% string for inline variable substitution

texts:
I'm using model builder to iterate through feature classes in datasets. With each sucessive process step (each model) I'm using the %Name% variable to name new feature classes in a new dataset.

I found that I can't just use %Name% without adding a little to the name of the new feature class; the model didn't seem to like working with two feature classes with the same name. To write the new feature classes to disk, I'm using something like %Name%_AppendedBit. The problem is that with sucessive steps I keep adding more and more bits to the name of the feature class.

I tried using Calculate Value with something like Replace(%Name%,_AppendedBit,_JustIdentifyThisStep) to make a variable with which to name the new feature class, but %output_value% doesn't appear to have any value (just "1"). I tried this as a variant and as a string. I was using vb Replace in the "expression" box, not messing with the "code block" box.

I could use %n% to name the feature classes, but then I lose useful detail from the original feature class name. This could be worked around by doing a spatial join later.

So, the long and the short of it is I'd like to use some kind of string truncation with inline variable substitution.


texts: Point data does not allign with polygon data, even though they are in the same projections

texts:
I created a feature class from XY data. This then adds a shapefile on my map, under the projection of WGS 1984. I then try to reproject my shapefile in to  Europe_Albers_Equal_Area_Conic. This works, but my shapefiles do not align with my already Europe_Albers_Equal_Area_Conic polygons. I click zoom to layer and they are far apart from each other. I have tried defining projection and reuploading the data, but nothing seems to work. I need all data in Europe_Albers_Equal_Area_Conic so that I can retrieve the correct depth measurements.

EDIT: Here is how far they are, even though they are projected the same, Figure 1 (picture).


texts: Incorporating Length of Route Value into Cost Evaluator Function (ArcGIS Network Analyst)?

texts:
A friend and I (who more frequently than I posts on this board) are creating a walking network and we have run in to a snag. We would like to attribute a cost evaluator to our network dataset using a self determined safety metric. We have the equation set up so that all we need to do is divide a field value based on the length of a route to get our "safety value". We just need to know in the cost field how to write the vb script or other evaluator so that Cost Value=[FieldValue]/"total length of individual route". Each route length will be different once each analysis is run, but we have a length value for each line segment, so we know there should be a way that adds up these pieces for the route for a total route length that can be divided by our field value. It seems like it should be a simple script line or something to put in the field evaluator box. Any suggestions?

Thanks


texts:
I have two tables with complementing data in them. Both have an ID column (tree_ID). Table A is missing several features which are in table B.

How can i copy only features from B to A which are NOT already in A using QGIS?

In pseudocode i search for something like:

Select features from B where B.ID is not already in A.ID


texts: "Undefined variable from import" Earth Engine and Python API

texts:
I'm writing some code for a GEE Application and I want to calculate the area of pixels in a given area.
I've found this method:

areas = (imNDVI4.multiply(ee.Image.pixelArea()).reduceRegion(reducer = ee.Reducer.sum(), geometry = test2, scale=10))

Unfortunately, I'm unable to run this part because of this error:

  Undefined variable from import: ...

that pops up both on pixelArea() and Reducer methods.

How could I solve this issue?

P.S. My "earthengine-api" package is already up to date.


texts: Creating web map with dynamically rendered symbols for multiple large datasets with ArcGIS?

texts:
I am editing an existing ArcGIS Server (v10.2) / Javascript (3x) map-based web app, which has several dynamic map services, some of which contain over 50,000 points.

The new web app must be able to render symbology based on a user-chosen attribute.  The user can then choose the map layer to view with a slider, choose the map service with buttons, and choose the attribute with another set of buttons.

Currently, I have a demo set up which uses a UniqueValueRenderer to create the renderer to pass to a LayerDrawingOptions object, which then is assigned to a given layer in a dynamic map service.

This is working somewhat, but I have:

some map network load errors,
the symbols rendered look kind of jagged and pixelated,
I haven't yet found a way to set scale levels for the symbols, besides tying every zoom event to re-running all the functions I am using to render the points (ScaleDependentRenderer won't work), and
I need to re-publish the services every month to account for new layers, which I do now with a Python script - but I haven't found a way to set the necessary properties to allow modification of dynamic map services which I do now in Server Manager manually for my demos (I posted define per request modification of layer order and symbology using ArcPy? about that concern)

Should I reconsider using feature services?

I was avoiding that because:

I didn't want to rewrite even more code and
from everything I've read, feature services can be sluggish when using services with many thousands of records.

If not that, is there any other way I can approach this problem?

At this time, I must use File GDBs and feature classes btw; so using a query based ArcSDE platform won't work.


texts: two symbols for a single label in legend

texts:
I am working on a spatial plan where I'm supposed to have two states - the existing one and the future one. I have always worked with separate maps - a single map for a single state. But I was wondering is it possible to put those two states together, i.e. to have two symbols in legend for one label - one on the left side of label and one symbol on the right side.
So far, I have one solution - to make two legends with appropriate symbols, put those legends close enough each other and turn off one of labels. But, is there a simpler way to do this?
I hope you understand my question.
Greets!

Here's an image of this problem.


texts: Listing the raster layers on the canvas for selection by user

texts:
I am creating a processing script that will take input from the user. The input parameter is to be a raster layer. Using the method described at https://www.qgistutorials.com/en/docs/3/processing_python_plugin.html
which lists vector layers that are on the canvas with
self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr('Input layer'),
                [QgsProcessing.TypeVectorAnyGeometry]
            )
        )

However, I want to have a source of raster layers. So I am looking to do exactly what the code does, but for raster layers. I cannot find how to do this. It seems from what I have found that it is necessary to list the filenames that the rasters are stored in. Is this the case? Can I list the raster layers for selection without  doing that?
I am using QGIS 3.16


texts: Is there a geoprocessing tool that I can use to find the closest road to a point on a map?

texts:
I have a Network Dataset layer that I use to perform routing.  Even if a starting or ending point is not near a road in the Network Dataset layer, I still want to perform routing using those points.
Is there a way that I can get the routing to start at the nearest place on the ND, to the starting point, or end on the route at the closest place near the ending point?


texts: Need a combination of geoprocessing tools that could replace Extract By Mask for use with WPF Runtime

texts:
I am in need of a combination of geoprocessing tools that could replace Extract By Mask for use with WPF Runtime.  The "Extract by Mask" tool is one of the tools that is not supported with Runtime.  Also "Reclassify" is one tool that is not being supported by Spatial Analyst extension.

However, right now I just want to concentrate on the Extract By Mask.  When I use that tool in my model, run the model and then try to create a geoprocessing package (gpk) from that result, I get an error message letting me know that that tool is not supported.  Therefore I need a combination of tools that could be used in my model, that would accomplish the same goal.

Thanks,
Renee


texts: GeoServer (CQL) query to return just the boundary polygon of a filtered feature

texts:
I have a polygon layer (in GeoServer). The goal is to just get the bounding box of a filtered feature. From the (E)CQL docs I thought the boundary filter should do the work, but I could not get it to work. Other threads found on Google didn't even help...

http://localhost:8080/geoserver/wfs/?service=WFS&version=1.1.0&request=GetFeature&typename=topp:states&cql_filter=BOUNDARY(the_geom)&

I always got the following exception:

java.lang.ClassCastException: org.geotools.filter.function.FilterFunction_boundary cannot be cast to org.opengis.filter.Filter
org.geotools.filter.function.FilterFunction_boundary cannot be cast to org.opengis.filter.Filter


texts: QGIS Processing Modeller - Select date by expression

texts:
I have a data file that contains current meter data for multiple days.

I am trying to use the Processing Modeller to extract a single days worth of data for further manipulation, using the extract by expression algorithm.

My input data layer has the date in it, with the format dd/mm/yy (i.e 01/12/19)

I've tested my selection expression using Select by expression in QGIS as per below, which correctly selects 1 days worth of data.

Building my model, I select the input later using a vector feature parameter definition.  I then used String parameter definition to input the required date '25/12/19' and the date field "%Date"

I entered the below expression into the Processing Modeller to match the select by expression that I tested ealier.

When I run the model, nothing is extracted to the result layer.

I'm sure it's probably a syntax problem, but cant figure it out.


texts: arcpy.Clip_management throws 999999 error. Running from toolbox works fine

texts:
I want to clip multiple rasters to my study area, so I'm writing a script that's using the arcpy.Clip_management tool. However, the tool throws an 999999 error with no explanation when executing.

I suspect that I'm mucking up the syntax somehow since the tool runs fine when executed from the toolbox with the same inputs, but I can't figure out what I might be doing wrong.

Tool help for reference:
http://desktop.arcgis.com/en/arcmap/latest/tools/data-management-toolbox/clip.htm

This is the script as is:

import arcpy, numpy

arcpy.env.overwriteOutput = True

try:
# set workspace
arcpy.env.workspace = r"C:\Users\Johnnie\Desktop\exarb\gis\pilotstudie"

# set inputs
inLocation = "in.gdb/"
outLocation = "ut2.gdb/"
studyArea = inLocation + "studieomrade"
vectorInputs = ["by_get", "my_get", "vl_get"]
rasterInputs = ["hojd2m3006_6589746_685960"]

# process: clip all vector inputs to study area
#for i in vectorInputs:
    #outName = i + "_clip"
    #arcpy.Clip_analysis(inLocation + i, inLocation + studyArea, outLocation + outName)

# process: clip all raster inputs to study area
for i in rasterInputs:
    outName = outLocation + i + "_clip"
    inName = inLocation + i
    arcpy.Clip_management(inName, "", outName, studyArea, "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")

except Exception as err:
print err.args[0]

This works with no problems for the commented out arcpy.Clip_analysis part.

This is what it looks like when I try to run the tool directly from the python console:

>>> arcpy.Clip_management("in.gdb/hojd2m3006_6589746_685960", "", "ut2.gdb/test2", "in.gdb/studieomrade", "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>   File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 13594, in Clip     raise e ExecuteError: ERROR 999999: Error executing function. Failed to execute (Clip).

I've tried with both empty quotes and "#" for default values as in the help sections code examples, but seems to make no difference.

Running the tool from the toolbox with the same inputs works just fine.


texts: Is there is a signal in QGIS when a project is about to be closed?

texts:
Is there is a signal in QGIS emitted when a project is about to be closed (QGIS  closing or some other project opening)?

I want to programmatically save changes in layers if they are in edit mode.


texts:
I have a multiline shapefile and a point shapefile. The points are supposed to fall in the center of line breaks but when imported to qgis, the points fall on one end of a line break. Is there a method to center these points? There is way to many to offset individually or selecting several points and using 'move feature' tool.
Here is what QGIS looks like and what the CAD drawing looks like. Im on QGIS 3.10.14


texts: Get info from WMS layer with OpenLayers mobile jQuery?

texts:
How to do a GetFeatureInfo request on a WMS layer (GeoServer)?

I am using OpenLayers jQuery mobile. I've already called the WMS in jQuery mobile, but how to get feature information when clicking the attribute.


texts: Merging various road data into single layer using ArcGIS for Desktop?

texts:
My company has two different road layers:

1) County roads
2) Old CAD roads

Each with their own pros and cons. I'm attempting to merge these into a single road layer that will be maintained.

The County has pretty good road data as far as location. It is also M-enabled for routing, which is nice. But it's missing road names.

Our CAD roads are dated and haven't been updated in a long time, but all the roads do have names.

My plan is to use the County roads as the base layer and try to join the CAD road names to it. But I first need to get them to line up correctly.

I believe the Integrate tool will work, but I haven't used it yet.

But then I need join the CAD road names to the County roads and keep the output M-enabled.

Any tips?


texts:
I have a SpatialPolygonsDataFrame and I would like to create an adjacency matrix. However, I would like to only count two polygons as adjacent if they share more than just vertices (i.e. they share an edge, not just a point). gTouches says polygons are adjacent if they share at least one point (e.g. just one vertex).

My thought was to run:

AdjMat <- gTouches(map, byid=TRUE)

then go look at the intersections for the adjacent polygons that identifies to see if the intersection is a line or a point.

test <- gIntersection(map[x,],map[y,])
exists("test@lines")

Playing around I find "test" is "Formal class SpatialLines" if x and y share an edge, and "Formal class SpatialPoints" if they just share a vertex (and "NULL (empty)" if they aren't adjacent at all. The problem I'm having is that even when it is "Formal class SpatialLines" and "test@lines" is actually in my environment, exists("test@lines") returns FALSE. (When it is "Formal class SpatialPoints, "test@lines" is not available, but "test@coords" is, but exists("test@coords") also returns FALSE.)

How do I get an adjacency matrix that doesn't count polygons as adjacent if their shared boundary/intersection is only one-dimensional?


texts: Which method to find intersection between lines with polygon

texts:
How could you find intersection between lines intersected with polygon?

I want to find the lines only that located inside the map and intersected with polygon then change the colour of intersected lines to red colour and extract the id of each line intersected in Excel table.

Edit: I used polygons to lines tool but it does not show me line and polygon layer, it shows me only one layer called Muffe and Arrow  as u can see in picture.

My file is dxf, I have imported by import/convert DXF then I got 3 layers text, line and polygon and Arrow as picture.

My question at first: Does it correct to import line and polygon at the same layer or not, because I though it will give me both of them seperately, one layer for line and one layer for polygon.

Then, when I used polygon to lines. So, I think, it considers polygon is line also so I do not need to convert it, am I right.

Then I used intersection tool and I choose line and polygon input and also for overlayer, for sure I got intersection will be the same input. And I am not sure what should I do.

Could someone guide me, how can I find intersection between lines inside the map and the polygon?

Something, I want to know, when I imort any file should I have seperately layers I mean, one layer for line, one for polygon and so on or both of them will be together.


texts: Adding Double field doesn't maintain precision in ArcMap

texts:
I have a shapefile that has an ID field(type=double, precision=15, scale=0). I am trying to add a new field(ID_2) that has the same field attributes as ID(Double,15,0). However, when I add the ID_2 field it seems like the new ID_2 field is being rounded to a precision=12.

I've tried adding the new field both manually within ArcMap Desktop 10.5 and also using the Add Field tool in a Python script; both methods result in the ID_2 field being rounded and the last 3 digits replaced with zero's. I also imported the same streets shapefile into a file geodatabase, added the ID_2 field, and calculated from the ID field, and that worked exactly as I would expect.

Also strange is if I right-click on the new ID_2 field and go to Properties, it shows the Precision=0.

Does anyone have any idea why my Double field isn't being added with a precision of 15, and is instead seemingly being rounded, chopping off the last 3 digits of the ID?


texts:
This is my first assignment for a GIS class, so bear with me. I have voting data for the state of California, and one thing I need to do is join a precinct voting data layer with a counties layer. When I click Validate Join, I get the following error message:

"Field names that match reserved words should not be used in database schema and can cause the join to fail. The following fields match reserved words;
- [TYPE] from < PrecinctVotingData >"

I know it's a simple fix but I'm stuck at this point. What should I do?


texts: Error updating Shape.area and Shape.length

texts:
am building a webservice with ArcObjects in C# .Net to insert a new feature or update an existing feature to a SDE-table. I am using the IFeatureCursor and an IFeatureBuffer to insert the features. When the geometry field ("Shape") is set, an COMException is thrown saying:

"Unable to update area or length field in response to shape update"

Obviously this refers to the default feature class fields Shape.area and Shape.length that are calculated automatically. I have also tried to insert features using the CreateFeature() and feature.Store() method, but this resulted in the same exception.

Just to be clear, this happens before I call the cursor.InsertFeature(featureBuffer) or cursor.UpdateFeature(feature) method.

My code looks like this:

private void Create(Result result)
    {
        var featureClassLoad = (IFeatureClassLoad) _master;
        var schemaLock = (ISchemaLock) _master;
        IFeatureCursor cursor = null;
        IFeatureBuffer featureBuffer = null;

        try
        {
            schemaLock.ChangeSchemaLock(esriSchemaLock.esriExclusiveSchemaLock);
            featureClassLoad.LoadOnlyMode = true;
            cursor = _master.Insert(false);
            featureBuffer = _master.CreateFeatureBuffer();

            foreach (var attr in result.Attributes)
            {
                var attributeName = attr.Field.Name;
                var index = _master.FindField(attributeName);
                if (attributeName == _master.ShapeFieldName)
                {
                    featureBuffer.Value[index] = attr.NewGeometry;//Exception thrown
                    continue;
                }

                if (attributeName == "ID")
                {
                    featureBuffer.Value[index] = "{" + result.Guid.ToString().ToUpper() + "}";
                    continue;
                }
                featureBuffer.Value[index] = attr.NewValue;
            }
            cursor.InsertFeature(featureBuffer);
            cursor.Flush();
        }
        finally
        {
            if (featureBuffer != null) Marshal.ReleaseComObject(featureBuffer);
            if (cursor != null) Marshal.ReleaseComObject(cursor);
            featureClassLoad.LoadOnlyMode = false;
            schemaLock.ChangeSchemaLock(esriSchemaLock.esriSharedSchemaLock);
        }
    }

Any suggestions are more than welcome!


texts: Snapping polygon layer to another layer automatically

texts:
I have two polygon layers in QGIS 3.4: 'Layer A' and 'Layer B'. These two layers have been digitized separately. Some features in these two layers are crossing or not meeting each other in some areas which I have to fix this problem.
'Layer A' is an accurate layer in terms of boundary so I want all the features in 'Layer B'  be snapped to the closest feature (point) in 'Layer A'.
Currently I try to move and snap all the vertexes in 'Layer B' to 'Layer A' manually using digitizing tool in QGIS which is very time consuming and I might miss some of them
So my question is: How can I do this process automatically or at least semi automatic?


texts:
I am trying to add mosaic/merge image bands with different layers on QGIS 2.18.0, however when I click the option to add mosaic I get the error:
ERROR 4: `' does not exist in the file system,
and is not recognized as a supported dataset name.

The entire stack trace specifically is:
File "C:\Qgis\bin\gdal_merge.py", line 540, in
sys.exit(main())
File "C:\Qgis\bin\gdal_merge.py", line 422, in main
ulx = file_infos[0].ulx
IndexError: list index out of range

Some questions imply it maybe empty spaces in the path string but, I manually chose the directory that has no spaces and even then I still get the error. How can I fix it?


texts: Uniform distribution of dots on world map depending on latitude

texts:
I would like to plot a map looking something like this:

I'd like to then connect each dot with the corresponding country and use it to display statistics. Therefore I need a set of coordinates, ideally calculated in python based on actual coordinates.

However, I want the dots to be uniformly distributed and I don't really know how to achieve that.
What kind of projection could I use? Is this even possible?

I tried a basic algorithm with a lot of "if" statements, which adjusts distance between two points based on the latitude, but the result is far from uniform. Basic code example:

def generate_lat(resolution):
    lat = 0
    lats = [lat]

    while lat < 90:
        if abs(lat) >= 80:
            step = resolution / 5
        elif abs(lat) >= 70:
            step = resolution / 4
        elif abs(lat) >= 50:
            step = resolution / 2.5
        elif abs(lat) >= 45:
            step = resolution / 1.5
        elif abs(lat) >= 30:
            step = resolution / 1.3
        else:
            step = resolution

        lat = lat + step if lat + step < 90 else 90
        lats.extend([lat, -lat])

    return sorted(lats)

def generate_coordinates(resolution):
    longitude_range = arange(-180, 180 + resolution, resolution)
    points = []
    latitudes = generate_lat(resolution)

    for lat_ind, lat in enumerate(latitudes):
        for lon in longitude_range:
            # offset every other line
            if lat_ind % 2:
                lon = lon + resolution / 2
            points.append((lat, lon))
save_file(points, resolution)


texts: Performing weighted sum overlay analysis in QGIS?

texts:
I have weighted raster data converted from polygons with different categories.

I have ArcGIS with 3D Analyst extension, unfortunately do not have Spatial Analyst, and I would like to use QGIS to make weighted sum overlay analysis, probably using ordered weighted averaging.

I found some similar questions but nearly all explained vector overlay analysis. Here grids (tif, img etc) will be used.


texts: Exporting all Map Series pages to PDF with page name as output file name using ArcPy

texts:
I am trying to export all of my current map series pages into multiple PDFs with the Page Name as the file name. So, "PAGENAME.pdf"... Currently, the code runs, but no output is given. I have looked around at other documentation and it doesn't seem to help in my instance. The output is not producing any PDFs at the moment but there are map series records in the project!
   import arcpy, os, sys

    aprx = arcpy.mp.ArcGISProject(r"\My Documents\20-0181 - Map Edits.aprx")
    outpath = (r"\My Documents\Phase 7 Data")
    l = p.listLayouts()[0]
    if not l.mapSeries is None:
      ms = l.mapSeries
      if ms.enabled:
        for pageNum in range(1, ms.pageCount + 1):
          ms.currentPageNumber = pageNum
          print("Exporting {0}".format(ms.pageRow.CCTVLines_Select_FACILITYID))
          pageName = ms.pageRow.CCTVLines_Select_FACILITYID
          l.exportToPDF(outpath + "{ms.pageRow.STATE_NAME}.pdf")


texts: Error with Intersecting Rasters and Lines in R

texts:
I'm trying to ''cut'' a spatial lines data frame object by a raster grid using R,

raster::intersect(SLDF, grid)

However, when I run this code I get,

Error in as.vector(x) : no method for coercing this S4 class to a vector

Here is my session info,

> sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.4 LTS

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
 [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
 [1] bindrcpp_0.2.2    viridis_0.5.1     viridisLite_0.3.0 flexsurv_1.1
 [5] survival_2.42-3   raster_2.6-7      rgeos_0.3-27      rgdal_1.3-1
 [9] sp_1.3-1          magrittr_1.5      dplyr_0.7.5

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17       pillar_1.2.3       compiler_3.4.4     RColorBrewer_1.1-2
 [5] plyr_1.8.4         bindr_0.1.1        tibble_1.4.2       gtable_0.2.0
 [9] lattice_0.20-35    pkgconfig_2.0.1    rlang_0.2.1        Matrix_1.2-14
[13] mvtnorm_1.0-8      gridExtra_2.3      grid_3.4.4         tidyselect_0.2.4
[17] deSolve_1.21       mstate_0.2.11      glue_1.2.0         R6_2.2.2
[21] purrr_0.2.5        ggplot2_2.2.1      scales_0.5.0       splines_3.4.4
[25] assertthat_0.2.0   colorspace_1.3-2   quadprog_1.5-5     muhaz_1.2.6
[29] lazyeval_0.2.1     munsell_0.5.0

Another strange thing, if I type raster::intersect into the console I get

> raster::intersect
standardGeneric for "intersect" defined from package "base"

function (x, y)
standardGeneric("intersect")
<environment: 0x5cc78f8>
Methods may be defined for arguments: x, y
Use  showMethods("intersect")  for currently available ones.


texts: GeoServer polygon transparency (SLD style file)

texts:
I am trying to add transparency to a set of polygons (city bounds) within a .sld style file for OSM.  Please see below:

As you can see the city bounds are overlapping Montana State University.  My goal is to have no fill or transparent fill for the city bounds. In this image I have it set to #00FFFFFF and GeoServer is not liking it as you can see.

Here is the SLD chunk for the polygon:

<sld:Name>Boundary</sld:Name>
<sld:Title>Boundary</sld:Title>
<sld:PolygonSymbolizer>
    <sld:Fill>
        <sld:CssParameter name="fill">#00FFFFFF</sld:CssParameter>
    </sld:Fill>
    <sld:Stroke>
        <sld:CssParameter name="stroke">#6a7c87</sld:CssParameter>
        <sld:CssParameter name="stroke-width">0.5</sld:CssParameter>
        <sld:CssParameter name="stroke-linejoin">round</sld:CssParameter>
    </sld:Stroke>


texts:
*Python 3.7 and QGIS 3.6

I have a program that loads multiple layers, and now I need to select different layers to edit. Before I was managing fine with iface.activeLayer(), but now that I have multiple layers this is no longer an option.

My current code - additional layers are loaded later:

layer = iface.addVectorLayer(in_file, name, 'ogr')

iface.mapCanvas().refreshAllLayers()
layer = iface.mapCanvas().layer(0)

layer.startEditing()

This results in the following error:

  AttributeError: 'NoneType' object has no attribute 'startEditing'

Interestingly, if I run the code when there is already a layer loaded, it works. This isn't an actual fix though, as it's not realistic for me to have a layer pre-loaded each time I run the code.


texts: Resampling multiple coarser raster images based on another finer resolution raster using R?

texts:
I have 2 raster images of different spatial resolution and projection.

I want to resample multiple coarser images (0.25 degrees) to finer resolution(4 km) based on one finer resolution image which is of 4 km. Both are in different folders. I am not getting how to do that.

I had written the code in R but I am getting an error.

library(resample)

>setwd("D:/mtech project/data/DEC 16/RH/1-DEC-16/extracted/")
>inFiles1 <- list.files(pattern="*.tif")
> nFiles1 <-  length(inFiles1)
> setwd("D:/mtech project/data/extracted/")
> inFiles2 <- raster("3DIMG_01DEC2016_0000_L2C_FOG_FOG.tif")
> for (i in 1:nFiles1) {
+     r1<-raster(inFiles1[[i]])
+     r2<-inFiles2
+     rs<-resample(r1,r2,method="bilinear")
+     write.Raster(rs, paste0('D:/mtech project/data/DEC 16/RH/RESAMPLED/rs_',i,'.img'), overwrite=T)
+   }

Error in .local(.Object, ...) :

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  :   Cannot create a RasterLayer object from this file. (file does not exist)

`


texts: PyQGIS: how to filter feature attributes according to symbol in the legend

texts:
I can reach feature attributes as follows:
for feat in layer.getFeatures():
    print(feat.attributes() )

This gives me whole attributes.
I want to filter attributes according to each symbol in the legend
Note that I can reach legend symbol items as follows
for symbolItem in layer.renderer().legendSymbolItems():

but I couldn't find any way to get attributes which represented by this symbolItem objects that have type of QgsLegendSymbolItem


texts:
I am using a road map in shapefile format in order to calculate a route between two points. The map is downloaded from OpenStreetMap  and it has small unconnected lines that are messing up my routing.

Is there any way that I could get rid of those without having to manually remove them one by one?
The lines look like this: http://screencast.com/t/4jnQFx57


texts:
I newbie have two newbie questions about Mapnik.

Does Mapnik allow generation tiles with float zoom level?
I have osm file for my region. I want generate images. Where i can download full xml file for it?


texts:
In QGIS 2.18.3, I want to set the alias for the fields.

layer.dataProvider().fields().field(0).setAlias("test")

That line returns no error but nothing is done.

I tried layer.updateFields(), but no alias set.


texts: Using zip function into distance calculating program?

texts:
I try to calculate distance between points from gpx file using harvesine formula. I stuck in moment when I should use 'lats' and 'lons' for formula which I found in internet. I want to use that two list as pairs, so zip function would be ok. But how i should do that?

import xml.etree.ElementTree as ET
import math
f = ('St_Louis_Zoo_sample.gpx')
p = ET.parse(f)
root = p.getroot()
lats=[]
lons=[]

for mainElement in root.findall('{http://www.topografix.com/GPX/1/1}wpt'):
        y = float(mainElement.attrib['lat'])
        x = float(mainElement.attrib['lon'])
        lats.append(y)
        lons.append(x)

def distance(origin, destination):
    lat1, lon1 = origin
    lat2, lon2 = destination
    radius = 6371 # km

    dlat = math.radians(lat2-lat1)
    dlon = math.radians(lon2-lon1)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c

    return d


texts: Mapbox GL JS - Trouble Loading Data to Create Layer

texts:
I am new to Mapbox GL JS and working with starter code from "Create a Time Slider" below:

https://docs.mapbox.com/mapbox-gl-js/example/timeline-animation/

To pull in the data, the code uses d3.json to perform a fetch (also I am new to d3....). I have imported this starter code into a Ruby on Rails app that I am creating for some basic data visualization.

https://docs.mapbox.com/mapbox-gl-js/assets/significant-earthquakes-2015.geojson

My goal is to modify the starter code to suit my own data set. So I had saved down the data...

https://docs.mapbox.com/mapbox-gl-js/assets/significant-earthquakes-2015.geojson

..into a local file earthquakes.geojson in the root directory of my app, assuming that the d3.json would just as easily point to my internal file instead of the Mapbox hosted URL

So the original starter code here:

d3.json('https://docs.mapbox.com/mapbox-gl-js/assets/significant-earthquakes-2015.geojson', function(err, data) {
if (err) throw err;

Is modified with a new URL for my local file:

d3.json('./earthquakes.geojson', function(err, data) {
if (err) throw err;

After making the above modification, none of the data pulls through, but I am not receiving any errors or info to point me in a new direction to investigate the issue.

Can someone give me a hint or point me in the right direction as to how to get Mapbox to connect with the data using d3?


texts:
I have a series of text files that represent X,Y,Z elevation data. I'm trying to use the GDAL grid command line tool to convert these text files to raster GeoTIFs. My files are space delimited and so I have created VRT files that point to the data so that GDAL can read them OK. The VRT files work. I can load the layers into QGIS and I can use ogrinfo to describe the data OK.

Many of my text data files convert to raster just fine. These are typically longer files with > 10,000 lines of data. However, my small files (< 5,000 lines) produce a raster with only NoData values.

I am using the linear algorithm. My points are evenly spaced at 1m and I've tried several radius values. The command that I'm using is:

C:/OSGeo4W64/bin/gdal_grid.exe -ot Float32 -of GTiff -txe 239982.0 240207.0 -tye 647659.0 648667.0 -outsize 225 1008 -a linear:radius=1 D:\Temp\19900928.vrt D:\Temp\19900928.tif

I've also tried the nearest neighbour algorithm and get some weird, partial values in the output GeoTiFF for these small files.

A sample source text file is available here:

https://gist.github.com/philipbaileynar/c9eb76880b8ce6fa8c786ca53d892b26

The corresponding VRT file is:

<?xml version="1.0" ?>
<OGRVRTDataSource>
    <OGRVRTLayer name="elevation_points">
        <SrcDataSource>CSV:D:\Temp\elevation_points.txt</SrcDataSource>
        <SrcLayer>elevation_points</SrcLayer>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:2762</LayerSRS>
        <GeometryField encoding="PointFromColumns" srs="EPSG:2762" x="field_2" y="field_3" z="field_4"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

Can anyone help and figure out how to make GDAL grid produce a raster from these data please?


texts:
I am working with Arcgis 10.3 and currently this problem keeps me awake at night. I import a csv file into ARCGIS with XY-coordinates and the Age of a person. I create a feature class which i later join with a layer file containing county boundaries. What i would like to know is how many Persons in a county fall in a certain age class.

My Goal is to get a table that groups the persons in a county in various age groups.

Name       Age0-30    Age31-60     AgeOver60
county1     25            10           9
county2     12            54           60

Is this possible and how?


texts:
I want to be able to draw labels with line. How can i do this?

Just like labels with line of some USA states in Google Maps, labels with red line in right :

Example from Google Maps


texts: ArcGIS javascript API: A way to position Dynamic Service layer on top of Graphic Layer?

texts:
I have 3 layers:

Tiled based map
DynamicServiceLayer with ~1000 points
Graphics layer with ~20 polygons

I really need the DynamicServiceLayer to be on top of the graphics layer (point must be on top of polygons), but I think it might not be possible - does the graphics layer always draw on top of all other layers?  I tried to reorderLayers without success.

Anyone drawn a dynamic service on top of a graphics layer?


texts: Batch Convert Folder of Shapefiles to KML / KMZ

texts:
I am looking to convert a folder with multiple shapefiles (over 100) to KML/KMZ files. I am sure I can build a model in ArcGIS to do this conversion but I am unsure on the iteration I should use.

I just need the KML files to have the same name and attributes as the shapefile. Style does not matter


texts: Add GPX files to postgreSQL database and add one column with a given value

texts:
I already found out how to add a gpx track information into my postgres database:

 `ogr2ogr -append -f PostgreSQL "PG:dbname=your_db user=xxxx password=yyyy" filename.gpx`

Everything works fine but now I'd like to have an additional column for each track, so that I can find a specific track...

Is this somehow possible?

A hack would be to execute a psql command after each import like this

`"psql -U postgres -d tracks -c 'UPDATE tracks SET track_id = 1234 WHERE ogc_fid = (SELECT max( ogc_fid ) FROM tracks);'"`

But maybe there is an option for ogr2ogr to do that?

Thanks


texts: Filter geoJSON layer after it has been loaded

texts:
I'm declaring and defining a geoJSON layer first, like:

var layer1 = L.geoJson(null, {
   onEachFeature: function(feature, layer) {
   ...
   }
});

..and then loading a CSV into it, like:

omnivore.csv('pbstatus.csv', null, layer1 )
   .on('ready', function(e) {
       ...
});

Now, I want to create other geoJSON layers as subsets of layer1 based on the value of one of the CSV's columns loaded into the layer's feature.properties dataset. Is there a way to filter layer1 AFTER it's been defined and loaded? I'm aware of the 'proper' way of filtering when defining the layer:

var layer1 = L.geoJson(null, {
   onEachFeature: function(feature, layer) {
   ...
   },
   filter: function(feature, layer) {
   return (feature.properties.foo == 'bar');
   }
});

.. but I don't want to repetitively do a different declaration each time. And, the data in the CSV is going to decide which layers and how many to form; we may not know all the possible filter conditions before we have imported the CSV. Is there a way I can spawn new layers as subsets of this layer1 after it's been loaded with geoJSON data ? Would the new layers also inherit the properties, actions and behaviours set in layer1's onEachFeature: function?


texts: How to access Esri ArcPy modules from Enthought Canopy?

texts:
I have ArcGIS 10.1 running on my windows system and  I have also installed the new 32-bit Canopy Express Enthought Python Distribution.

Is it possible to load and use the Esri ArcPy modules from Canopy?


texts: Layer Styling by Geometry Generator in QGIS 3.6

texts:
I am trying to create custom styling by using Geometry Generator in QGIS 3.6 but I got an error expression.

segments_to_lines ( $geometry)
 ramp_color('Viridis',
  scale_linear(
        length(
            transform(
     geometry_n($geometry ,
     @geometry_part_num),
     'EPSG: 4326' , 'EPSG: 26910'
     )
     ) / (
     m(end_point(
     geometry_n($geometry ,
     @geometry_part_num)))

    -
    m(start_point(
     geometry_n($geometry ,
     @geometry_part_num)))
     ) * 3.6,
     0,50,
     0,1
     ))


texts: Extraction of reflectance value in ERDAS Imagine

texts: How can I extract the reflectance value from a raster in ERDAS Imagine. If possible I would like to know the steps.


texts: I want to create a line style in QGIS (model is attached here with),


texts: Mapnik renders almost blank image when rendering 2 layers

texts:
From Python code I can get Mapnik to render this PNG image:

And this one:

But if I try to render both these layers (to render the dots on top of the basemap), then I get only this PNG - which is all transparent; but showing a "half dot" at the bottom. Almost is if the image is rendered but the clipping is wrong:

There are no error messages or other clues. What can I be doing wrong? I sense something fuzzy with unaligned projections; but I am not sure.

The Python code looks like this:

import mapnik
mapnik.logger.set_severity(mapnik.severity_type.Debug)

# Base map
map = mapnik.Map(1600, 1200)

# Create raster and point symbolizers
raster = mapnik.RasterSymbolizer()

point = mapnik.PointSymbolizer()
point.file = 'dot.png'
point.allow_overlap = True

# Create new rules and add the symbolizers
raster_rule = mapnik.Rule()
raster_rule.symbols.append(raster)

point_rule = mapnik.Rule()
point_rule.symbols.append(point)

# Create new styles and add the rules
raster_style = mapnik.Style()
raster_style.rules.append(raster_rule)

point_style = mapnik.Style()
point_style.rules.append(point_rule)

# Add the styles to the map
map.append_style('Raster style', raster_style)
map.append_style('Point style', point_style)

# Specify our data sources
raster_datasource = mapnik.Gdal(file='K100_2017_1_20_612_62.tif')
point_datasource = mapnik.Datasource(type='csv',file='stats.csv')

# Create new layers for the map, add the data sources and styles to those layers.
raster_layer = mapnik.Layer('Raster')
raster_layer.datasource = raster_datasource
raster_layer.srs = '+proj=longlat +datum=WGS84 +units=m +no_defs +over'
raster_layer.styles.append('Raster style')

point_layer = mapnik.Layer('Point')
point_layer.datasource = point_datasource
point_layer.srs = '+proj=longlat +datum=WGS84 +units=m +no_defs +over'
point_layer.styles.append('Point style')

# Add the layers to the map

# If both of these two lines are active - a blank image is rendered?
# If just one of them (any) - and image is rendered fine?
map.layers.append(raster_layer)
map.layers.append(point_layer)

map.srs = '+proj=longlat +datum=WGS84 +units=m +no_defs +over'
map.zoom_all()
mapnik.render_to_file(map, 'K100_2017_1_20_612_62-base-pos.png', 'png32:z=1')


texts: Python script in ArcMap Field Calculator giving Error 000539: SyntaxError: Invalid Syntax?

texts:
I've been trying to calculate this expression for a field in an attribute table in ArcMap 10.4.1, and keep getting the following error:

  Error 000539: SyntaxError: invalid syntax (expression, line 5)

60 *math.degrees((math.acos(math.sin( math.radians(!LAT_O! ))
*math.sin(  math.radians(!LAT_DEST!) )
+ math.cos( math.radians( !LAT_O!) )
*math.cos( math.radians( !LAT_DEST! ))
*math.cos( math.radians(!Dist_Long! )))

They are all number fields - Any idea where I'm going wrong?


texts:
I have GeoServer (v2.2 from OpenGeo Suite) running on a test server hosting many different WMS services sourcing many different PostGIS db tables.  I want to push this setup onto a production server now but I'm not sure what the best method would be for migrating this info without recreating all the Workspaces, Stores, and WMS layers again.

Both servers are the same:

Windows Server 2008 R2 Standard

My thoughts so far are:

Stop GeoServer on test
Make backup of all dbs in PGAdmin on test
Restore backup(s) to production db
Copy data directory from test and paste it to same location on production
Start GeoServer on production

Thoughts?


texts:
I would like to use a startup.py which is part of a rollout and not located inside the users config-folder (like ~/.qgis2).

Plugins located in the apps folder (e.g. C:\OSGeo4W64\apps\qgis-ltr\python) are working but if I try to use a startup.py in this folder it isn't executed.

Is there something I can do to get this startup.py executed or is this a bug?


texts:
Can anyone give me direction as to possible reasons why when publishing certain map services I get the "Packaging Succeeded but Publishing Failed" error?

When I expand the Geoprocessing Results window the failure is when trying to Upload/Publish the Service Definition file.  If I try to rerun the process from the .SD file I get an error message telling me that my .SD is not supported or doesn't exist; when I follow the path to the folder location my .SD is actually a .SDDRAFT.  If I change the file extension to .SD the process will run but it just fails like it did originally.

The frustrating part is that I have no problems publishing some things.  For example, a map service with only vector data publishes with no concerns, an image service with LiDAR in a GDB publishes without an issue, JPG rasters publish without issue how ever TIFF rasters are failing, any form of aerial imagery in a GDB is failing, TIFFs in a GDB are failing.  There seems to be no rhyme or reason to it as the same dataset will publish for one project but fail for another.


texts:
I want to calculate the elevation of a particular feature. I have both DSM and DTM as input data, both with the same cell size and same number of columns and rows (same resolution). For the DSM, min = 532.807 and max = 541.081. For the DTM, min = 532.771 and max = 537.959.

I want to calculate the elevation (DSM - DTM) of each pixel and finally add all of them to calculate the final elevation of the .tif file representing the particular feature. I followed the procedure given below using GDAL in QGIS.

dtm = gdal.Open(r"C:\Task\Data\DTM.tif")
dsm = gdal.Open(r"C:\Task\Data\DSM.tif")

# taking pixel values into a numpy array
dtm_nparray = np.array(dtm.GetRasterBand(1).ReadAsArray())
dsm_nparray = np.array(dsm.GetRasterBand(1).ReadAsArray())

# numpy array representing Height of each pixel
height_nparray = np.subtract(dsm_nparray, dtm_nparray)

# DSM - DTM = Elevation
sum_elevation = np.sum(height_nparray)

Unfortunately, I got the value of sum_elevation as a negative value. This means that DSM - DTM = -ve value which is actually wrong.

What is the exact issue here ?


texts: Symbology in QGIS: raster with complex discrete data

texts:
I have a raster with the following pixel values:

1100: corresponds to grassland
1104 to 1134: corresponds to a forest plantation and values from 4 to 34 are the age of the plantation
1204 to 1234: corresponds to another type of forest plantation and the values from 4 to 34 are the age of the plantation
1304 to 1334: corresponds to another type of forest plantation and the values from 4 to 34 are the age of the plantation
And many more values of this type.

How can I assign colors with the symbology to this raster like this?

1100: yellow
1104 to 1134: shades of green from low to high intensity
1204 to 1234: shades of orange from low to high intensity
1304 to 1334: shades of red from low to high intensity

It is not an option to transform the raster into a vector because it corresponds to a very large area.


texts: I am trying to figure out flow direction for a stream network and I do not have a high resolution DEM data for the same. I only have information about the stream orders and the network itself. Is there a way to figure out the flow direction based on stream-order connectivity as an inbuilt function?


texts: Weighted flow accumulation with pygeoprocessing

texts:
I cannot manage to get results for a weighted flow accumulation with pygeoproccessing. When I try it without the weight, it does work. However, I need the weight and when I include it, python just stays in calculation mode forever without useful results. It does create a raster file for the flow accumulation, but without accumulation results. Also, there is no error message in the python console. One time, I let it continue a whole night and nothing more came out.

I prepared the dem in QGIS with fill sinks (saga) and unified extent and resolution of the dem and the weight raster with QGIS raster calculator, so the dimensions of the datasets are the same.

I have tried it with data for a complete basin, but also with a small clipped test area. The same problem occurs.

My code (working in anaconda, spyder, python 3.7):

#load module
import pygeoprocessing.routing as geop

#define test data locations
path = "C:/test/"
dem = "{}dem_clip.tif".format(path)
flow_dir = "{}flodir.tif".format(path)
ro = "{}ro_clip.tif".format(path)
flow_acc =  "{}flowacc.tif".format(path)

# calculate flow directions
geop.flow_dir_d8((dem,1),flow_dir)

# accumulate
geop.flow_accumulation_d8((flow_dir,1), flow_acc, (ro,1))

I usually do the weighted flow accumulation with saga (in QGIS/R), which works, but since im running into issues with larger datasets I am trying to explore other ways. That is the reason I am trying pygeoprocessing.

Does anybody know what can be the issue?


texts: QGIS / "create grid" consumes far more memory on linux than on windows

texts:
I am using QGIS 3.10 both on windows and on Linux with the same RAM and CPU resources.
When running the tool "create grid" QGIS behaves strange on Linux: For a certain extent (~ 10x10km), when creating a grid with 20m-100m cell size, memory consuming doesn't go beyond ~2,5-3GB. With 15m cell size, the calculations seem to exponentiate a lot - and it takes up to 30GB memory (not a steady increase, but rather abrupt) and it won't finish either. As my goal was initially a 10m grid - with this cellcize QGIS seems to crash a few minutes after running the tool, as memory goes up abruptly, however I don't any messages or warnings.
On windows, running the same tool with the same parameters behaves differently: for the 10m grid only half of the memory (max. 7-8GB from 16GB) is used - however I didn't let it complete as the area still seems to be too large.
I am splitting up the area for the calculation of the grid now, but I still wonder what causes this strange behavior on Linux ?


texts: QGIS3 - Simple automation of QGIS output - turn values on/off & print out dozens of maps

texts:
I need to make a series of maps of a single location and make a single map for each value in the layer.
I have a layout I've created with borders, texts, north arrow, and everything, with the map, and an overview map on it.
My current workflow is:

Select value 1.
Go to print composer
Output to image
Go back to map
turn off value 1
turn on value 2
go to print composer
output to image...

So on and so forth for 34 images. It seems like something that would be easy to do with automation, but I have looked at tutorials for atlas, and the timecontrol export, but haven't seen anything that does this with layouts and that turns on and off values.
I have been able to make the atlas work with separate layers, but I don't know how to copy the color and symbology over so it match up from the values to the layers, without doing it 1 by 1, defeating the whole purpose of saving time & automation.
Can you give it a tutorial link or video or point me in the right direction?


texts: Read only a specific value in LC_Type1 Class Table for MCD12Q1.006 MODIS Land Cover Type Yearly Global 500m

texts:
I have defined a specific area to look at but now I want to only read the data for the water.

MCD12Q1.006 MODIS Land Cover Type Yearly Global 500m has a LC_Type1 Class Table and I only want to read value 17 because it has the water. How can I do this?

var Landcover = ee.ImageCollection("MODIS/006/MCD12Q1"); //terra
  var collectionTerra = Landcover
  .filterBounds(Lake_Sup)
  .filterDate('2014-07-01', '2014-08-31');

  var clip1 = function(images) {
    return images.clip(LakeSP);
  }


texts: Converting Formal Class RasterLayer to Large RasterLayer?

texts:
I should note I'm using RStudio. These terms may be specific to RStudio.

I've seen a couple of questions referencing this, but I haven't found any explanation of how they are different. I'm finding calculations are taking a very  long time with a formal class rasterlayer but not a large raterlayer, so I want to understand how they are different and how to solve this.

Here is an example raster. It's only about 2mb. When I load this raster, it loads as a formal class rasterlayer. However, if I load this raster and perform any simple mathematical operation to it, it becomes a large rasterlayer. For some reason, though (and this is not happening with this example raster), this takes about 10 seconds to perform any simple calculation with my data.

r1<- raster("./example_raster.tif")
r2<- r1 + 0  #this takes about 10 seconds with my data

Is there any other way to convert from the formal class rasterlayer to a large rasterlayer?


texts:
I am using QGIS and I simple want to clip a shape with another shape but the result isn't right - the result is only that what you see in my screenshot. What can I do?

Picture 1: Shape 1

Picture 2: Shape 2

Picture 3: Result


texts: WFS - from PostGIS to Shapefile with GeoServer

texts:
From server side, I'm using GeoServer and from client side I'm using OpenLayers to allow the users see the layers (PostGIS) on the map.
The new requirement is that the users may download layers in shapefile format.

What is the best way to do this?

I would prefer add the option to download them in the same map view.


texts:
Sorry for the lack of terminology as it escapes me.
I have a shapefile that contains 96000 individual polygons and I'd like to simplify the shapefile so that it creates 1 polygon.  I think it is something like the Region or Aggregate, but only for shapefiles.
Reason being, I have 100s of shapefiles and by using the merge command creates a huge shapefile that contains many smaller polygons of similar features, but I'd like a shapefile to contain just one large polygon.  Think of Hawaii being in the middle of the ocean, but it is part of the US, so when selecting any US state (Montana), it also selects the rest of the states along with Hawaii.  I also want to retain the original geometry of the individual polygons.
I am running the merge tool right now and it has been running for over two hours and my shapefile is at 2 gb and counting.


texts:
I'm looking into geocoding U.S. locations without using Google or Bing.
I found a service that geocodes according to the Tiger Census data of 2013. I have little knowledge about this Census and didn't get extremely far poking it around.
When I use the api of the service that geolocates based on that census, I get several possibilities for the same full address in some cases. E.g. I got that one address in San Francisco is related to 3 different zip codes (that address, randomly picked, is "1354 Market St, San Francisco, CA")

Is that at all possible in the U.S. zip code system?

Do you think the Tiger Census data of 2013, is a very good source for geocoding?

Thanks,
Matan


texts: Is there some standard way of obtaining bounding box of an S-57 chart, e.g., for one downloaded from https://www.charts.noaa.gov/InteractiveCatalog/nrnc.shtml?


texts: Vertically docking Tools toolbar to left of map and right of TOC in ArcMap 10.x?

texts:
In general, I really like the docking environment in ArcGIS 10.  As I was showing off ArcGIS 10 to some of my end-users, I couldn't get the toolbar docked to the left of the map and to the right of the table of contents (TOC). Am I missing something here?  I want to configure like ArcMap 9.3.1:

I think as a workaround, end users can dock the toolbar to the left of the TOC, and auto-hide the TOC(personally, I will do this).  Or they can dock the TOC to the right of the map (and maybe also auto hide the TOC).

To quote @MattWilkie:

  Please vote for ArcMap to bring back vertical toolbar docking, next to
  table of contents at
  http://ideas.arcgis.com/ideaView?id=087300000008GlOAAU


texts: Uploading DEM heightmap in QGIS with CRS that gives coordinates in degrees and correct distance measurements

texts: I have a DEM heightmap (1600px * 800px) showing a fictional planet in equirectangular projection that has the same dimensions as Earth. I want to upload it in QGIS with a CRS that gives coordinates in degrees and correct distance measurements. The best result I've managed to obtain is a map where the coordinates (shown in the right bottom corner of the QGIS window) go from -180°...+180° east-west, and +90°...-90° north-south like I wanted, but where the distances aren't correct (the "measure line" tool gives a horizontal distance of 40000km at the equator, which is correct, but also 40000km at the poles, where it should be 0km).


texts:
I'm using Django and Openlayers to add some data from my PostGIS database to the map but can't seem to work out how to style the vector layer that uses the kml.

var layer = new ol.layer.Vector({
    title: 'Polygon',
    source: new ol.source.Vector({
      projection:new ol.proj.get("EPSG:3857"),
      url:'http://localhost:8000/layer',
      format: new ol.format.KML(),
      extractStyles: false,
    }),
    });

    map.addLayer(layer);

I've been through lots of tutorials attempting to use ol.style.Style and styleMaps and more but whatever I do the layer stays as just white. I was wondering if anyone knew if this was just a recent change in the OpenLayers code meaning there aren't any tutorials or if there is a solution anyone may know about.


texts: Sentinel-2 remove images from image set with B10 above a value (remove cirrus)?

texts:
I've noted in Sentinel-2, some of my images in a collection appear hazy (and have B10 values > 50) yet the cloud coverage assessment shows no cloud/cirrus.
I understand I can see the cloud_coverage_assessment for the whole image, what what about for a specific geometry (ROI)? What if we wanted to limit the images in an imageset to cloud free conditions in the ROI only (could be completely cloudy outside the ROI, I just want cloud free within the ROI)?
(If one clips an imageset, i.e.
var images = imageset.map(function(image) {
return image.clip(ROI)
})
and wanted to remove images from the imageset that had suspected cloud cover (B10 values above some threshold) within the ROI, how to go about this?


texts: Can I make an ESRI checkout replica from a checkout replica?

texts:
If I create a checkout from a clients geodatabase, can I then turn around and put that replica on our ArcGIS Server and create a couple replicas from it?

So, then, our replicas would be synced to the "replica" on our server which would then be exported back out to a personal geodatabase.

Would that personal geodatabase be able to be synced back to the clients database or would doing all this ruin what makes it a replica?

I hope this made sense.

Basically, what I'm trying to accomplish is getting a replica from a client that multiple people will be working on, but before syncing it back to the client, we would like the replicas we worked on synced back to 1 database so we can check the connectivity before giving it back to the client.


texts: Is there any WCS service example with Mapserver?

texts: I want to publish tiff image in browser using Mapserver, for obtaining elevation data. Is there any resource or example regarding this.


texts: Some data loss when Importing Personal Geodatabase Table from Excel?

texts: When import table from excel some data were loss. I was trying to change field type but I can't solve this error. I have attached an example image to show the problem.


texts: ImportError: cannot import scipy from windows powershell

texts:
I am having an importError message while trying to import scipy from powershell. But I am able to do this when I try to do this in the python shell in ArcGIS.
The error I get is this:

Traceback (most recent call last):
File stdin, line 1, in module
ImportError: No module named scipy

While trying to import arcpy, numpy, gdal all works well. Can anyone help with the problem? I did search but could not find a similar problem here.


texts: Options for displaying PostGIS rasters in OpenLayers

texts:
This question is similar to an earlier question about displaying vector data in a web-map. I want to have a basic web interface which uses OpenLayers and can display raster data that is currently in PostGIS, using a basemap such as google. I have installed PostGIS2.0 with the raster library and see from the documentation that there are several options to output data, such as JPEG, GeoTIFF and PNG. Is their a recommend format to use?

To display vector data I am currently using javascript to send coordinates from the user's browser to my server, and then performing a query using PHP/SQL on PostGIS. The results are then returned as GeoJSON which can be overlaid on the basemap. I wondered if a similar approach could be applied using rasters but I don't know what raster format is best, or if there are any other constraints I should consider.

I have seen a lot of references to GeoServer but I'm not sure if it is necessary for this project (I'm not entirely sure when GeoServer should be used or not).


texts: Rotate symbols perpendicular to line angle (ArcGIS)

texts:
I have lines representing curving streams (blue line in the Picture) and want to angle symbols perpendicular to the angle in their particular Locations (point features symbolized as red lines in the picture).

With Layer Properties - Symbology -> Advanced I can use a field Expression to rotate my Symbole as desired.
I therefore Need to create a field in the attribute table with the correct angle-number in it which I then can sue in the field expression. Does anyone know, how I can formulate an Expression which tells me the angle the Point has to rotate to be perpendicular to the line?

I work with ArcGIS 10.2.

Desired:


texts: Converting River Channel Cross-Section Data to Bathymetry

texts:
I have multiple surveyed cross-sections of river channels and need to convert this data into 3D polylines (shapefile) for use in a hydraulic model. The cross-section data is typically supplied in "X,Y,Z" csv format.
I have tried using this method described by Purdue University, but it fails to run each time except when I use their input example datasets. I also have limited access to ArcMap. I've tried a method described by modelling software "HEC-RAS", but this requires significant data conversion to formats compatible with HEC-RAS. From my experience, it also produces an output that does not follow the river banks correctly and often contains interpolation errors.
I plan to write a python script that will project the first cross-section onto a second and interpolate intermediate points within each cross-section so that both have elevations at corresponding points. Between each corresponding point, I will then use the river centre line and river bank lines to inform interpolated polylines along the channel connecting the corresponding points in each cross-section. I should end up with something similar to the images below, except curved to match the river banks and centreline.
Is there a simpler way to interpolate this data into 3D polylines along the river channel? Ideally, any methods will use QGIS 3+ or just python outside of a GIS, but I do have access to ArcMap 10.4, MapInfo 10.5 and Global Mapper.


texts: What does the vector to raster option in QGIS do, exactly?

texts:
I have a large point vector layer (a few million objects) with measured elevation values and have converted it to raster, using the raster>conversion>(vector to raster) option in QGIS. It seems to work fine, but I haven't been able to find info on how it exactly the value for a pixel is computed.

Does QGIS take the average of all points in the cell, the median, the first object it finds in the attribute table, etc.?

Is there perhaps a plugin with options for this?


texts: How can I export a Seed Points layer as a image?

texts:
I have the Seed Points layer produced by Watershed Segmentation. I want to export is as a image (with background).

How can I do this?

I understand that I can do Project->Export->Export Map to Image.

However this takes the extent of the points layers, which might be smaller than the image to which the points have been put. When it should be the extent of the original image.


texts: Map scale will only double and halve in QGIS

texts: When I zoom in or out in QGIS, my scale will only halve or double respectively. This is the case whether I manually type the scale into the window at the bottom of the screen or select one of the set scale options. The scale automatically jumps to the nearest scale that is an exponent of base 2,256 (the smallest scale allowable for some reason), where the exponent is a multiple of 2. For example, if my scale is 144,447 and I type in or select 250,000, it automatically jumps to 288,895. This is occurring prior to adding any layers. I probably just have to uninstall/reinstall but though I would check here first.


texts:
I am developing a GIS application and in my code I need to call the exiftool executable which is installed in a different path.
I call in command prompt and in OSGEO4W shell and both are fine.
But when I write in my code doesn't give any error but doesn't work.

import os
import subprocess
os.chdir('C:/micmac/binaire-aux')
comand3 = 'exiftool.exe (...) args (...)'
subprocess.call(comand3)

In the command lines I use this code:

C:/micmac/binaire-aux/exiftool (...) args (...)

The args are exactly the same.

When I try just:

import os
os.chdir('C:/micmac/binaire-aux')
comand3 = 'exiftool.exe'
subprocess.call(comand3)

This works fine. I don't understand why doesn't work with the arguments, if they are the same.


texts:
I am reading xml data with FME 2014 and attempting to rebuild a relational database model. I have an xml object "element " which has an m:m relationship with object "type". The xml looks like this...

<element id="194070" >
  <types>
    <type id="445"/>
    <type id="643"/>
  </types>
</element>

Using fme I can flatten the xml and recreate a table which looks like this...

element.id | type.id.1 | type.id.2
194070          445         643

But in order to build a association table for the m:m relationship I need this;

element.id | type.id
194070          445
194070          643

Can anyone offer a way to create this association table using fme transformers?


texts: Which Geo stack for an online strategic wargame?

texts:
I'm googling around for some weeks but I've some difficulties to figure out where to begin my project and which part of the available Geo stacks are mandatory for my needs.

I have a strong Java background, some interest for Python and totally open to learn some new technologies.

As an hobby and because it is close related to my job, I plan to develop a web MMO strategic wargame heavily based on an interactive maps (basically to display the tactical situation and give orders to units).

Here is what I was thinking about :

Geoserver to serve maps and layers.
PostGIS to store geodatas.
GeoExt to display in a browser
Apache Click as the Java application server (combat resolution, unit management)

I managed to deploy a GeoServer (war in tomcat, easier is difficult)
I managed to install PostGreSQL
I run some tutorials about GeoExt. I need to experience Javascript a little bit.
I played with Apache Click to serve pages and forms.

Here is the flow I was thinking of :

GeoServer serves the map backgrounds and layers constructed from PostGIS data for static objects (Cities, Victory Points, Territories, ...)
Apache Click generates HTML+GeoExt pages in the client browser.
Apache Click sends CRUD commands to PostGIS in response of user's action.
[something] updates regularly PostGIS (movements, combat resolution, unit characteristics management...)
Apache Click generates dynamically KML files to live update the browser (KML NetworkLink ?) for units display with icons, zones (range, detection...)...

And now, the missing brick in my little wall :

How can I implement a "sea/land" layer. For example, to invalidate a ship movement passing on land ? GeoServer layer ? Zones defined in PostGIS ?
Is the following flow reactive enough to update the map after a user action : HTML/GeoExt event -> Apache Click process -> PostGIS update -> GeoServer layer update -> HTML/GeoExt display

The main point is : where are the data if I want to use it for calculation and update ? GeoServer static layer or PostGIS tables ?