texts: Seeking open source tools to visualize map projection distortion?

texts:
A reference article on Map Projections contains the following figure:

I think the images are created using Geocart, listed in Additional Resources section. I'm interested in trying to create visualizations like this--particularly the areal distortion and angular distortion visualizations--using open source tools. Either desktop GIS or code (Python, R, etc.) would be helpful.


texts: Adding continuous color ramp as legend using ArcMap

texts:
I am interested in adding a continuous color ramp legend with numerical labels like the following example to place on the layout view of arcmap:

I am quite new to arcMap so I'm not sure how I can do this. On the symbology tab, I used graduated colors with a preset color ramp from ArcMap. There are 60 classes, so when I insert a color ramp, the result becomes this:


texts: Creating Isochrones - pgr_drivingDistance 2.6 SQL error

texts:
Im trying to create a matrix of driving distances from an osm2pgrouting parsed dataset.

Following many blogs I have the below incorrect SQL:

SELECT * from ways join (
SELECT
*
FROM
pgr_drivingDistance(
    'SELECT gid as id
        , source
        , target
        , length_m as cost FROM "ways"
        '
    ,274738,500,FALSE)
) as route on ways.gid = route.id

I get the error  " ROUTE.ID" does not exist.  This seems sensible since the pgRouting docs for 2.62 show no field returned in the examples called ID.

If I replace route.ID with say route.seq (a column PG driving distance does return) i get a very dodgy result with points random scattered over the map.

What is the equivalent of ID or what do i need to change to get a proper matric returned of points with 500m?


texts: Unable to find inherited method for function ‘crop’ for signature ‘"numeric" using R

texts:
I'm trying to iterate a crop and then a resampling process by using a for loop but I've faced some problems because I'd obtain an error message. Below is my code, then the error message
#find all tifs in your directory
dir<-"E:\\Proyecto PM2.5\\2_PM_2.5_Processing\\Test"

#get a list of all files with .tif in the name in the directory
files<-list.files(pattern='.tif', full.names = TRUE)

# raster with the expected extent
d <- raster("E:\\Proyecto PM2.5\\3_PM_2.5_Entrega\\temporal\\Raster_D.tif")

#raster with the expected characteristics: extension, cellsize, number of pixels
r_ref <- raster("E:\\Proyecto PM2.5\\3_PM_2.5_Entrega\\temporal\\Raster_C.tif",crs="+init=epsg:32614")

for (i in 1:length(files)){
  r <- raster(files[i])
  rq <- crop(r[i],d)
  rp <- raster::resample(rq[i],r_ref,
                 method = "bilinear",
                 filename = file.path ("./res", files[i]),
                 format = "raster",
                 overwrite = TRUE)
}

By running this I obtain the following error:
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘crop’ for signature ‘"numeric"’

I've read that could be related to that I'm applying an S4 method to a not supported object, but looking into the class for r I obtain this:
[1] "RasterLayer"
attr(,"package")
[1] "raster"

How do I solve this error?


texts: When to enable 'On the fly' CRS transformation in QGIS? When to turn it off?

texts:
I understand what 'On the fly' does, and I have re-read https://docs.qgis.org/2.6/en/docs/training_manual/vector_analysis/reproject_transform.html

I tend to have it turned on all the time, but I don't know if that is correct or not? I am aware it can cause issues.

I have had a look through the related questions and they all tend to be when users are having problems, and further complicated by people asking about their particular datasets and respective CRS's.

Is there simple rule-of-thumb answer to when to have 'On the fly' turned on and when should it be turned off?


texts: Field Calculator not showing some fields

texts: When I open my attribute table, click on a column, and select Field Calculator it opens and doesn't list all of the Fields in my attribute table. Why can't Field Calculator find these Fields?


texts:
I've created  a hillshade raster (ArcGIS 9.3) and I have 16-Bit tiff image (please see attached).
I need to submit 1 band/8-Bit with grey level (0-255) and set NoData for all pixels outside of the survey limit. My problem is while I'm changing to 8-bits, my NoData pixel changes to 0 value, which is wrong. How can I change this correctly and don't change NoData to 0?

For simple operation I can also use ArcGIS 10, I don't have any extensions for this version.


texts: Installing PostGIS behind firewall on server with no Internet connection?

texts:
I'm trying to install PostGIS on a server so that I can run some SQL scripts for styling in QGIS, but to do so it seems to have to install through Stack Builder which requires an Internet connection.

I have installed PostgreSQL, so how do I install the PostGIS plugin locally?


texts: Trouble Detecting State When Coordinates Near Coastline

texts:
Using Python I'm attempting to determine the state a coordinate pair resides within.

I've downloaded the US Census State And Equivalent dataset, extracted the shapefiles, and put them into Shapely objects.

sf = shapefile.Reader('shapes/tl_2017_us_state.shp')

for record in sf.shapeRecords():
    shapely_shape = shape(record.shape)

    # Ive also tried to parse the shapefile like this according to:
    # https://gis.stackexchange.com/questions/113799/how-to-read-a-shapefile-in-python?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa
    # shapely_shape = shape(record.shape.__geo_interface__)

    if isinstance(shapely_shape, Polygon):
        shapely_shape = MultiPolygon([shapely_shape])

    point = Point(37.907478, -122.369055)
    if shapely_shape.touches(point):
        print(record.record[5])

According to Google Maps the point (37.907478, -122.369055) should be on the California coast, however the .touches() doesn't appear to show the intersection.  My suspicion is that I'm losing some accuracy when going from a shapefile to a shapely object, so I'm looking for some insight into what I'm doing wrong.


texts: postgis values within XYZ viewport in openlayers

texts:
I got this XYZ server where I get my maptiles.

Now, when the user is at zoom level 1 ( closest). I want it to execute a command that gets the current viewport from the map and Queries my PostGIS database to get the geometry points that are contained within that viewport.

How could this get done? and what is the best way to display this data as fast as possible on the map?

Thank you in advance :)


texts: How do I align symbols of side by side columns in qgis 2.14 print composer

texts: Using QGIS 2.14. I created this legend in the print composer using the columns tab, but as you can see the top symbols don't align horizontally at the top. What is the setting to tweak to have these symbols align correctly? I feel like i've been through everything and can't find it.


texts: Using Field Calculator to create origin and destination lines on map?

texts:
I'm a novice GIS user.

I'm trying to create origin and destination lines on a map. I now I will have to eventually use the points to lines to tool. My question revolves around the attribute table and field calculator. I have to latitudes and longitudes of each location. I need to create a new field to show the transition between the origin and destination points. I am confused about how to go about this and what equation I should put in the field calculator.


texts: Copying reduceRegion() results to featureCollection properties

texts:
I'm calculating a time series for a 7-year period over a grid tessellation AOI. My goal is to have a CSV of zonal statistics for all features in a feature collection from 3 different bands in an image collection.
To do this, I'm running reduceRegion() on an image collection using a combined reducer, over a feature collection of 55,000 grid cells. The image collection has 3 different image bands (NDBI, NDVI, UI). I want to copy the six resulting reduction values (2 for each of the 3 bands) to each feature as properties, then export those values as columns in a csv. However, I am having trouble actually copying the results of the reduction to the feature collection as properties.
My Script (https://code.earthengine.google.com/b58961f43a176ea520400c36a4940de1)
What I have tried:

A nested .map() function that reducesRegion by each image in a collection, and then maps those results to the feature collection (line 125-136):

// combine the max, mean reducers
var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.count(),
  sharedInputs: true
});

var scale = ee.Image(l8sr_bands.first()).projection().nominalScale();

// try copying image bands as a feature property?
// this works w/o error, but no bands are copied to feature properties. reducer values must be returned or copied back to featureCol, but not sure where.
var l8_data = table.map(function(feature) {
  return l8sr_bands.map(function(image) {
    return ee.Feature(feature)
      image.reduceRegion({
        reducer: reducers,
        geometry: feature.geometry(),
        bestEffort: true,
        scale: scale
      }).copyProperties(image).copyProperties(feature);
  });
}).flatten();
print("output test 1: ", l8_data.first());

Alternate code that copies the six reducer results (mean and sum for all 3 bands) as properties to the images (line 149-162). This does not do exactly what I want, because it produces results as properties of an image collection. I want to map these properties to the feature collection.

var reducers = ee.Reducer.mean().combine({
  reducer2: ee.Reducer.count(),
  sharedInputs: true
});

var scale = ee.Image(l8sr_bands.first()).projection().nominalScale();

var getStats = (function(img) {
            var image = ee.Image(img);
            var reducers = ee.Reducer.mean().combine({
                reducer2: ee.Reducer.count(),
                sharedInputs: true
                });
            var stats = image.reduceRegion({
                reducer: reducers,
                geometry: table.geometry(),
                bestEffort: true,
                scale: scale
                });
            return image.setMulti(stats);
})

// map all stats onto the feature collection while copying the image properties as feature properties
// result = sum and mean for each band as image properties. how do I map this to the feature collection?
var stats_fc = table.map(function(feature) {
  return l8sr_bands.map(getStats);
}).flatten();
print("output test 2: ", stats_fc.first());

I believe I'm missing a return or feature.set() statement somewhere to copy the reducer values back to the features, but I'm not sure where to do this. I'm new to GEE and JS, so I'm sure it's something simple.
I am aware that there are similar questions already on GIS SE, such as this question and ReduceRegions not writting to feature properties in GEE, as well as the developer's forum, but posted solutions have not worked for me.


texts: Exporting render of polygon after applying smooth geometry command using QGIS?

texts:
I have recently discovered the smooth($geometry,x) command in QGIS which is very convenient for me. However, when exporting polygons, the "smoothness" obviously goes away as the vertices are untouched by the command. (the command is shown here: https://www.youtube.com/watch?v=yEurXCS1hEw )

Is there a way for me to extract the rendering of the polygon as it shows after applying the command?

The goal being to be able to load the polygon later on in other projects without having to reuse the command, and with vertices that match the render I had previously achieved.

I'm attaching a picture to help explain what I'm looking for.


texts: Erase Overlapping QGIS Polygon using Geopandas within Jupyter

texts:
I have a Polygon of transects that contain overlapping features where they converge. I need to erase where they overlap, but still keep one of the polygons there (preferably the preceding one).
Essentially, I want to iterate through each feature, have the current feature (f1) see where it overlaps with all other features (for loop, f2) erase that overlap from f2 and assign its new geometry as that.
I've used this code within the QGIS Python Console which works, but I'm developing a notebook to automate the generation of transcripts, and this is the final task. Here is the code I use within Q:
for f1 in layer.getFeatures():
    for f2 in layer2.getFeatures():
        if f1.id() < f2.id():
            geom1 = f1.geometry()
            geom2 = f2.geometry()
            new_geom = geom2.difference(geom1)
            layer2.dataProvider().changeGeometryValues({f2.id(): new_geom})

I found a chunk of code which gets me a shapefile of every overlapping area, now I need to erase that geometry, but only from the succeeding feature (f2) not f1. Here is the code.
g1 = gpd.GeoDataFrame.from_file(linesBufferFolder + "\\Lines_Buffer")
g1.shape

import itertools
geoms = g1['geometry'].tolist()
intersection_iter = gpd.GeoDataFrame(gpd.GeoSeries([poly[0].intersection(poly[1]) for poly in itertools.combinations(geoms, 2) if poly[0].intersects(poly[1])]), columns=['geometry'])
intersection_iter.to_file(linesBufferFolder + "\\Lines_int_iter.shp")

I'm assuming I need something in that for loop to tell it to erase while it is there?

Above is the end result. The highlighted polygon comes after the one overlapping it on its right, therefore the polygon to the right remains the same, and the highlighted on is now the shape you see. If I highlight the one to the left of the currently highlighted, it would also show no overlap. All of these polygons are square 1m buffers.


texts: Writing Numpy array to raster file (tif) returns a trivial black square

texts:
I need to turn a 2D numpy array into a raster file that I can later import in ArcGIS (version 10.3.1, licence type: Advanced). The coordinate system that I need to target is the British National Grid coordinates (OSGB1936, EPSG:27700), which are expressed in meters.

This is how I create the tif. This procedure is based on the answer by Eddy The B.

from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import numpy as np
import matplotlib.pylab as plt

#These 3 arrays come from a netCDF file
array=inv_masked
lat=lats #degrees
lon=lons #degrees

xmin,ymin,xmax,ymax = [130859.699132,-964.660232,655859.699132,527035.339768] #meters (obtained with a conversion)
nrows,ncols = np.shape(array)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)

geotransform=(xmin,xres,0,ymax,0, -yres)

output_raster = gdal.GetDriverByName('GTiff').Create('C:\Users\MyName\Temp\inv_masked.tif',ncols, nrows, 1 ,gdal.GDT_Float32)
output_raster.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(27700) #British National Grid OSGB1936
output_raster.SetProjection(srs.ExportToWkt())

The inv_masked array contains boolean values only, and has 352 rows and 350 columns. When plotted with matplotlib, this is what inv_masked looks like:

However, the .tif file that I create with such procedure is just a black square. The spatial reference seems to be correct, but when loaded in ArcGIS I don't see the coastline. Essentially, the .tif map only has 0 values, so all the 1 are not preserved, somehow, and that's why the coastline is not visible. What am I doing wrong here?

What am I doing wrong here?


texts:
What is the easiest way to view a raster in CartoDB?

From the list of formats that you can upload into CartoDB, I noticed that .tif is available, so I've tried adding a GeoTIFF, but it says that there is no geometry and it doesn't display anything.


texts:
I'm preparing a map with more than 50 sites listed. Such sites are labeled with numbers and I would like to add a text box at the bottom of the map with the full site name's list. However, it seems that there is no "justify text" option in QGIS Print Composer.

Does anyone has a suggestion on how to achieve this?


texts: Compressing multiple static layers to single layer to make web map loading faster?

texts:
In our web application there are many layers to be loaded, this cause the loading of map very slow.

I am seeking a way to compress multiple static layers to a single layer to make the loading faster.

All the layers from an enterprise geodatabase (version 10).

Is it possible to compress layers?


texts: Problem with units in distance matrix after projection change in QGIS

texts:
In QGIS 2.16.3, I have a point file, originally in the WSG84 (EPSG:4326) projection. I wish to determine the distance in km between each point in this point layer and intend to use the distance matrix tool to do this.
As WSG84 uses the unit dd, I first transformed my projection to what I believe is a metric transformation (Plate Carree (this was chosen as it was described in the QGIS manual as being proper when the goal is to measure distances)).

Problematically, when creating the distance matrix, the unit still comes out as dd, not km. I have also tried using the UTM transformation, with no luck.
After reading online, I saved my project (save as) and reopened it. Now, the drop down menu under Vector has changed, and I can no longer choose the distance matrix at all.

How do I calculate the distance between points in km?

I have now used the Save as function to save my point file in the UTM projection. But now the Vector drop down menu is altered, and I can no longer choose the distance matrix:


texts:
I'm looping through all the layers in my map using

map.getLayers().forEach(function(layer, i) {
   ...
});

I'd like to get the name of each layer but I haven't been able to find anything in the examples or the documentation. Any pointers?

I'm using OpenLayers 3.0.0.


texts: Exporting Graphic with transparency in layout to PDF using ArcGIS for Desktop?

texts:
I'm using ArcGIS 10.1 for Desktop.

I am trying to export a pretty simple map to PDF and finding the results are not what I expect. I have inserted two graphics on top in Layout. One is just a rectangle that has a custom user style applied. The style uses a Picture Fill derived from a PNG file that has transparency set to 80% opacity in the original file.

The other graphic is a rectangle text box that has a similar Picture Fill Style applied, this one using a 20% grey fill and 80% opacity in the original PNG file. On screen, the boxes look fine. That is, I can see other layer information underneath that looks greyed out with a partially transparent box overlaying.

When I export to PDF, the boxes show up as empty of any fill at all. The sample text in the text box shows up, but not the transparency fill.

Resolution is set to 300 DPI, Resampling is set to best. Colour space is set to CMYK. There is no compression. Bitmap fills are set to vectorize. I have tried varying every setting that can be changed in Export to PDF and find that nothing really changes with regard to this transparency issue.

Any ideas? I have scanned other forum posts and see that it is an issue from several versions back. Has it ever been solved?


texts: What kind of resampling will give the highest quality when building pyramids?

texts: The resampling setting of the pyramids builder of QGis povides several ways to do so : nearest, average, gauss, cubic, average_mp, average_magphase and mode. I am not a specialist in the subject. I just wonder what way of resampling renders the highest output quality, no matter the output size.


texts: Transform the whole ImageCollection into a normalized histogram for a region in Google Earth Engine

texts:
I'd like to transform the whole ImageCollection into normalized histograms for a region. Each image should have one histogram. I want to merge all the histograms into one Image (each column corresponds to a histogram of the original image) or a Feature, and then export the Image. It seems ee.Reducer.fixedHistogram can help, and I can do it for one image. But when I use the map function applies it to the whole collection, it doesn't' work.
var dataset = ee.FeatureCollection('TIGER/2016/Counties');

var region = dataset.filterMetadata('STATEFP', 'equals', '01');
var region = ee.FeatureCollection(region).filterMetadata('COUNTYFP', 'equals', '047');
var region = ee.Feature(region.first());

var imgcoll = ee.ImageCollection('MODIS/006/MOD09A1')
      .filterDate('2002-12-31','2012-12-31').select([0,1,2,3,4,5,6]);

var imgcoll2=ee.Image(imgcoll.first());

var countyImg = imgcoll2.reduceRegion({
  geometry: region.geometry(),
  reducer: ee.Reducer.fixedHistogram(-100,9000,32),
  scale: 500 // the resolution of the GRIDMET dataset
});
print(countyImg)


texts:
I have a python code where I am trying to assign a unique, random value to every row. I have tried:

maxvalue = int(arcpy.GetCount_management(shape).getOutput(0)) + 10000
for row in arcpy.SearchCursor(shape, fields="FID"):
     insert = randint(0, maxvalue)
     arcpy.CalculateField_management(shape, "RNDM", '"' + str(insert) + '"', "VB")

but to no avail. I keep getting only 1 value for all rows. I have looked at other resources, but the answer still eludes me (e.g. https://geonet.esri.com/thread/96917). This method is as well extremely slow.


texts:
I have been given some images in TIFF format, and I need to convert them to the PNG format. I have been given that the TIFF file contains 4 bands: R, G, B, and NIR (which I do not need for the PNG). Further, the TIFF file contains the pixels in uint16. On running gdalinfo -stats on a TIFF file, this is what I get:

Driver: GTiff/GeoTIFF
Files: 1.tif
   1.tif.aux.xml
Size is 1364, 1295
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
    AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
    AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",9],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32632"]]
Origin = (471456.599999999685679,5252603.999999998137355)
Pixel Size = (0.600000000000145,-0.600000000000000)
Metadata:
  AREA_OR_POINT=Area
  Band_1=PC Sharpen (Band 1:02AUG16102514-M2AS-005695178010_01_P001.TIF)
  Band_2=PC Sharpen (Band 2:02AUG16102514-M2AS-005695178010_01_P001.TIF)
  Band_3=PC Sharpen (Band 3:02AUG16102514-M2AS-005695178010_01_P001.TIF)
  Band_4=PC Sharpen (Band 4:02AUG16102514-M2AS-005695178010_01_P001.TIF)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  471456.600, 5252604.000) (  8d37'17.52"E, 47d25'34.53"N)
Lower Left  (  471456.600, 5251827.000) (  8d37'17.70"E, 47d25' 9.37"N)
Upper Right (  472275.000, 5252604.000) (  8d37'56.59"E, 47d25'34.66"N)
Lower Right (  472275.000, 5251827.000) (  8d37'56.76"E, 47d25' 9.49"N)
Center      (  471865.800, 5252215.500) (  8d37'37.15"E, 47d25'22.01"N)
Band 1 Block=1364x1 Type=UInt16, ColorInterp=Gray
  Description = PC Sharpen (Band 1:02AUG16102514-M2AS-005695178010_01_P001.TIF)
  Min=10.000 Max=1118.000
  Minimum=0.000, Maximum=1652.000, Mean=217.443, StdDev=89.146
  Metadata:
STATISTICS_MAXIMUM=1652
STATISTICS_MEAN=217.44318436577
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=89.146238735783
Band 2 Block=1364x1 Type=UInt16, ColorInterp=Undefined
  Description = PC Sharpen (Band 2:02AUG16102514-M2AS-005695178010_01_P001.TIF)
  Min=30.000 Max=1705.000
  Minimum=0.000, Maximum=2133.000, Mean=309.619, StdDev=145.897
  Metadata:
STATISTICS_MAXIMUM=2133
STATISTICS_MEAN=309.61859565892
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=145.89744142331
Band 3 Block=1364x1 Type=UInt16, ColorInterp=Undefined
  Description = PC Sharpen (Band 3:02AUG16102514-M2AS-005695178010_01_P001.TIF)
  Min=0.000 Max=1594.000
  Minimum=0.000, Maximum=2082.000, Mean=201.293, StdDev=126.679
  Metadata:
STATISTICS_MAXIMUM=2082
STATISTICS_MEAN=201.29347252573
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=126.67918628134
Band 4 Block=1364x1 Type=UInt16, ColorInterp=Undefined
  Description = PC Sharpen (Band 4:02AUG16102514-M2AS-005695178010_01_P001.TIF)
  Minimum=0.000, Maximum=3102.000, Mean=453.964, StdDev=208.556
  Metadata:
STATISTICS_MAXIMUM=3102
STATISTICS_MEAN=453.96443347411
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=208.55563519719

When I convert it via gdal_translate, I get a black image without the -scale option, as the maximum values are far from the uint16 maximum (65535). With the -scale option (and after removing the alpha channel using ImageMagick), this is what I get:

However, this is not the actual output, as what I can view in QGIS is this (image obtained by exporting the map; which does not keep the original tif image size):

How is it that QGIS gets the correct output? I have tried obtaining the numpy array format of the tif file and manipulating it to get the output to match what I get in QGIS in vain. Since I need to convert multiple TIFF files to PNG, I need to write a script to process all of the TIFF files. How do I do this?


texts: How do I sum up polygons in each squared grid in QGIS?

texts:
So I have created a square grid (fishnets) layers with another layer of Philippine mangrove distribution polygons, I intersected two layers so the Area1 (philipine mangrove polygons) and Area2 (square grids) and so their both area columns are in the same attribute table. However I want to calculate the percentage of mangrove polygons in each grid cells, I would need to sum the mangroves polygons for each grid cell.
How can I do that? I know a way of doing it is via ArcMap's Tabulate Intersection function. Is there another function equivalent to qgis?
I also tried with join attribute by location from the data management tools, however the output layer's attribute table would just not contain the sum of Area1. Is it because that the mangrove polygons are not each wholly contained in grid cells but intersecting? What should I do then?


texts: Change Map layers when zooming in (JavaScript API)

texts:
I have created a web map application (Javascript API) :

require([
  "esri/map",
  "esri/dijit/InfoWindowLite",
  "esri/InfoTemplate",
  "esri/layers/FeatureLayer",
  "dojo/dom-construct",
  "dojo/domReady!"],
function(
  Map,
  InfoWindowLite,
  InfoTemplate,
  FeatureLayer,
  domConstruct
){
  map = new Map("map", {
    basemap: "topo",
    center: [-85.416, 49.781],
    zoom: 5,
    logo: false});

  TestLayer = new esri.layers.ArcGISDynamicMapServiceLayer(
    "http://Myserver/arcgis/rest/services/TestMAP/MapServer");

  query = new esri.tasks.Query();
  query.returnGeometry = true;
  query.outFields = ["CCSNAME", "DAUID20111" ,"DAUID201_5" ,"DAUID201_9" , "DAUID20_14",  "DAUID20_19", "NAME" , "CITY",  "POST_CODE", "HRUID2007.csv.GPF", "HRUID2007.csv.GPE", "HRUID2007.csv.SPECF", "HRUID2007.csv.SPECE" ];

  symbol = new esri.symbol.SimpleFillSymbol(
    esri.symbol.SimpleFillSymbol.STYLE_SOLID, new esri.symbol.SimpleLineSymbol(
      esri.symbol.SimpleLineSymbol.STYLE_DASHDOT, new dojo.Color([0,255,0]), 2),new dojo.Color([0,255,0,0]));

  map.addLayer(TestLayer);

In the default zoom level set, 3 map layers from **http://Myserver/arcgis/rest/services/TestMAP/MapServer** are showing up.

How can I change my code so that in the default zoom level I get another 3 layers from this url:

"http://Myserver/arcgis/rest/services/TestMAP2/MapServer"

and then when I zoom in I get the current ones.


texts: WMS tile layer getfeatureinfo OpenLayers

texts:
I'm using this link to try and add the WMS GetFeatureInfo Tile Layer functionality to my OpenLayers map
http://astuntechnology.github.io/osgis-ol3-leaflet/ol3/05-WMS-INFO.html
My code is currently as follows:

// Create a popup overlay which will be used to display feature info
var popup = new ol.Overlay.Popup();
map.addOverlay(popup);

// Add an event handler for the map "singleclick" event
map.on('singleclick', function(evt) {

    // Hide existing popup and reset it's offset
    popup.hide();
    popup.setOffset([0, 0]);

    // Attempt to find a marker from the planningAppsLayer
    var feature = map.forEachFeatureAtPixel(evt.pixel, function(feature, layer) {
        return feature;
    });

    if (feature) {

        var coord = feature.getGeometry().getCoordinates();
        var props = feature.getProperties();
        var info = "
" + props.casereference + "
";
            info += "
" + props.locationtext + "
";
            info += "
Status: " + props.status + " " + props.statusdesc + "
";
        // Offset the popup so it points at the middle of the marker not the tip
        popup.setOffset([0, -22]);
        popup.show(coord, info);

    } else {

        var url = lakes_wms.getSource().getFeatureInfoUrl(
                        evt.coordinate,
                        map.getView().getResolution(),
                        map.getView().getProjection(),
                        {
                            'INFO_FORMAT': 'application/json',
                            'propertyName': 'NAME,AREA_CODE,DESCRIPTIO'
                        }
                    );

        reqwest({
            url: url,
            type: 'json',
        }).then(function (data) {
            var feature = data.features[0];
            var props = feature.properties;
            var info = "
" + props.NAME + "

" + props.DESCRIPTIO + "
";
            popup.show(evt.coordinate, info);
        });

    }

});

However, I'm getting the error expmap.js:565 Uncaught ReferenceError: request is not defined I've tried going back on the linked example, but they don't seem to define it either. So I'm not sure what I'm defining it as in order to continue


texts: Generating random number of points in each polygon based on polygon area?

texts:
I would like to generate random points in various polygons from the same layer in ArcGIS 10.2 (I tried the "create random points" feature), however I do not want the same number of points in every polygon.

I would like to set it so that the number of points per polygon is based on the size (with a minimum of at least 1 point per polygon), so that the bigger polygons have more points generated in them and the smaller ones have less.

I am doing this to figure out points to ground truth underwater benthic reef features.


texts: Pyramid Creation Utility for ArcGIS 10.1 or later

texts: Someone wrote a really neat ENVI Pymramid Creation Utility I am wondering if there is similiar one for ArcGIS


texts:
I assume the answer to the number of CRS is "nobody knows", given that everybody can define their own CRS willy-nilly.

So maybe it's better to ask: How many EPSG-registered CRS are there? I haven't found an easy way to scrape the total number of EPSG codes in the various lookup tools.

Are there other well-endowed libraries of CRS such as EPSG? If so, how many CRS do they entail?


texts: how to display polygons as points until certain zoom level?

texts:
can somebody tell me how to display polygons as points until you zoom to some level. I'm speaking about OpenLayers.
So, for example, if the zoom level is less then 10 than it should display feature as a point and when you zoom to 10 or further it should display it as polygon with the shape it really has.

I've tried something like this:

var my_style = new OpenLayers.Style({
'fillColor': '#9900cc',
'fillOpacity': .8,
'strokeColor': '#330066',
'strokeWidth': 1,
'pointRadius': 4
});
var vector_style_map = new OpenLayers.StyleMap({
'default': my_style
});

map.events.register("zoomend", map, function () {
    if (map.getZoom() < 10) {
        console.log("zoom level is less then 10");
        layer.styleMap = vector_style_map;
    }
});

But it just paint my polygon with defined colors, it doesn't change the shape to point.
I'm probably missing something obvious but I would appreciate if you could help me.
Thanks


texts: How to add a VFR map coordinate grid in QGIS?

texts: I want to add a coordinate grid like in the link below to my map. A simple grid is no problem, but I would like to know how to add the fragmentation of the lines.


texts: Is it possible to implement sticky move tolerance in a ArcGIS JavaScript API application?

texts:
We are working on a project in ArcGIS JavaScript and is going fine. We want to add sticky move tolerance to our project so it will not be so easy to move features on the map.

Is it not possible to add that to our project? Does anyone know anything about that?


texts: What are common methods for converting raster digital elevation data to vector contours?

texts:
In researching this I have found many references to the opposite, creating elevation data from contours, but nothing regarding contour creation. Preferably, I am looking for open source solutions like QGIS, Grass or command line tools but I am open to any suggestions that will get me pointed in the right direction.

I am planning to take the contour data into a vector drawing app for line formatting but if the app can do some of that as well it would be a major time saver.


texts: Currently using the Maximum Likelihood tool in ArcGIS Pro to pick up green leaf pixels in an aerially captured RGB image. Does the tool use an algorithm which filters pixels based on thresholds in the HSB (hue, saturation & brightness) color space? At present I can't seem to find any documentation supporting the algorithm breakdown.


texts: Time Enabled Rest Services & the Pro Time Ribbon

texts:
When consuming a time enabled REST service in ArcGIS Pro 2.7, the time extent does not seem to update or read correctly when using the the Time Ribbon. I am consuming the https://nowcoast.noaa.gov/arcgis/rest/services/nowcoast/radar_meteo_imagery_nexrad_time/MapServer service. This should have radar data for the past 4 hours. When I set the Full Extent to be that layer (as I understand it, the start and end times of that layer should then populate the start and end times for the full extent), it is inputting an incorrect date (for me it keeps setting it to some times on 1/26/2021).
This has happened regardless of whether I use the service directly via Add Path, or as a server added to the project or favorited. Refreshing the server had no effect. The same time window is put in for other time enabled services from that server.
The only solution has been to manually set the Start & End times of the full extent range, but this is not a useful work-around for such continuously updated data.
Does anyone have any idea about what is going on here? My only thought is that Pro cannot read the time extent of a REST service layer.


texts: Strip Map Index Data driven pages not creating in order

texts:
I have a map book I am trying to create from a line that i pieced together from 2 different feature classes.  I combined the 2 different lines, and made sure my end vertex is at the end of the line, but my map book is still placing the pages out of order.  I have tried creating the line several different ways, and was wondering if there is something I am missing.  I have tried it about 5 different times and it is not working.  I have created a strip map index off of this line before, and have not run into this problem.

To reference, the line is only one segment, and the end vertex is at the end of the line, there are no segments that cross over each other.


texts: Setting coordinate systems in gvSIG: workflow or transformation issues

texts:
I am learning how to use gvSIG 1.10 (Mac OS 10.6.8 Snow Leopard). My goal is to create a simple county map of Illinois using the base data shapefiles available from the Illinois Natural Resources Geospatial Data Clearinghouse.

I want to use a state plane coordinate system such as NAD83/Illinois (West) so that the layer will look like this (narrow):

... not this (wide):

Here are the dialogs I see after adding the layer:

I tried to change the projection to the desired "narrow" version, but no matter what I do, it is always displayed as the undesired "wide" version.

In case this is relevant, here's the empty Transformations list in gvSIG:

What am I doing wrong? Am I missing a software component?


texts: Using ArcGIS for Server and ArcSDE for importing shapefiles?

texts:
Is there a way to use ArcGIS Server web services to perform the task below?

I have .NET ArcObjects code to import shapefiles to Spatial DB which works on ArcGIS for Desktop.

This code used ArcSDE to parse several folders and import into Spatial DB.

What kind of ArcGIS for Server web service can I use to web enable this functionality?

Browser will have a login credentials and just needs to show progress of importing and failure/error conditions in browser.


texts: Joining a queried table to a feature class

texts: I have a table in a geodatabase that I would like to join to a feature class. The table contains well locations with multiple sample dates and associated data, however, I am only interested in sample events that occurred in August. I applied a definition query to the table, but when I tried to join it to the feature class it failed. I removed the definition query and tried to join again and it was successful, but joined the data was based on the earliest sample event (February). I am working with ArcGIS Desktop v10.3. Is it possible to apply a definition query to a table and join it to a feature class?


texts: Finding combinations of polygons in PostGIS

texts:
all,

I have found something that has made me a bit curious about how to correctly use ST_INTERSECTION(A,B).  I am looking for every combination of two political district boundary sets in California, the state senate and the state assembly. These are two different sets of districts in California, both of which cover the entire state.  There are fewer polygons in the upper chamber (senate) than in the lower chamber (assembly), which means that each upper polygon will be bigger and intersect with multiple lower polygons, possibly containing them.

The shapefiles that I am using are from the US Census, but I have preserved the specific state of the files that I am using in data.world if you want to try them.

The first thing to note is that loading these into PSQL via shp2pgsql results in them coming in as MULTIPOLYGON, so I recast them to POLYGON using these queries:

DROP TABLE IF EXISTS state_leg_lower_single;
CREATE TABLE state_leg_lower_single AS
  (
    SELECT
      sll.name,
      sll.statefp,
      sll.sldlst,
      sll.lsad,
      sll.lsy,
      (ST_DUMP(sll.geom)).geom :: geometry(Polygon, 4326) AS geom
    FROM
      state_legislative_lower sll
  );

DROP TABLE IF EXISTS state_leg_upper_single;
CREATE TABLE state_leg_upper_single AS
  (
    SELECT
      slu.name,
      slu.statefp,
      slu.sldust,
      slu.lsad,
      slu.lsy,
      (ST_DUMP(slu.geom)).geom :: geometry(Polygon, 4326) AS geom
    FROM
      state_legislative_upper slu
  );

Both of these seem to come through correctly in QGIS.

Lower Chamber Districts

Upper Chamber Districts

Now, here is what Quantum GIS does when I use Vector -> Geoprocessing -> Intersection, which is what I expected.

To get this result in PostGIS, I have tried so many variations on the following query:

DROP TABLE IF EXISTS state_legislative_intersections;
CREATE TABLE state_legislative_intersections AS
  (

    SELECT
      'SLL-' || sll.name || ': SLU-' || slu.name AS composing_districts,
      sll.name                                   AS state_leg_lower_name,
      slu.name                                   AS state_leg_upper_name,
      COALESCE(
          ST_INTERSECTION(sll.geom, slu.geom),
          ST_INTERSECTION(slu.geom, sll.geom)
      )                                          as intersecting_geom
    FROM state_leg_lower_single sll, state_leg_upper_single slu
    WHERE
      (ST_INTERSECTS(sll.geom, slu.geom) OR ST_INTERSECTS(slu.geom, sll.geom))
      AND
      GeometryType(st_intersection(sll.geom, slu.geom)) IN ('POLYGON', 'MULTIPOLYGON')

  );

What I get is a subset of what I imagine should be the result, only 71 polygons, with massive lacunae.

This has become a matter of curiosity to me, so I would like to see what I am doing incorrectly. Can anyone help?


texts: Updating only one of two maps in QGIS Print-Composer without updating other map

texts:
I have two maps within the map composer in qgis. Whenever I want to update one map it updates both maps.

How can I prevent QGIS from updating the other map?


texts:
I have a custom made Toolbox and I can't make it run.

The Toolbox is called toolbox_k and the tool is called sn

I did this:

arcpy.ImportToolbox(r"C:\Users\user\Desktop\VST\Toolbox_k.tbx", "sn")

I use the sn with the parameters as using the GUI version and it doesn't work. It gives a syntax error:

sn('pathone','pathtwo')

UPDATE

Running:

>>> arcpy.ImportToolbox("Toolbox_PS.tbx", "tbxk")

gives:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\__init__.py", line 125, in ImportToolbox
    return import_toolbox(input_file, module_name)
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\toolbox_code.py", line 457, in import_toolbox
    mymodule = generate_toolbox_module(toolbox, None, False, False, False, module_name, use_alt_alias)
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\toolbox_code.py", line 432, in generate_toolbox_module
    'exec')
  File "C:\Users\user\Desktop\VST\Toolbox_k.tbx", line 19
    def 00CalculateKAEK(PST=None):
         ^
SyntaxError: invalid syntax


texts: How to Access/Write to a Shapefile in S3 Bucket using BOTO3 and GDAL

texts:
I'm trying to access and make edits to a shapefile that's located in a zip on an S3 bucket. So far I've used /vsizip/vsis3/ to access the zip file that the shapefile in question came in. I've made the necessary edits in a MEM file in OGR, and now need to upload to S3, either in a zip or directly to the bucket.
What's the best way to accomplish this task? Sadly vsis3 does not allow for random writes. No preference on how, just need to get this data back to S3 bucket.
Edit: I'm using python


texts:
I would like to host an Openstreetmap server using Geoserver (Geoserver/GeoWebCache/postgis). However, I encounter a problem regarding the tiles and gridsets. TMS urls (myserv/z/x/y.png) are different than urls from openstreetmap.org or MapQuest. Furthermore, there is an offset on the y axis.
The following url shows an area from my own server. The other url shows the same area from Mapquest.

http://devmobicite.dnsroute.fr:8080/geoserver/gwc/service/tms/1.0.0/mapServer:mapOSM/14/16280/12392.png

http://mtile01.mqcdn.com/tiles/1.0.0/vy/map/15/16280/11637.png

The offset seems to be caused by the zoom level (here 14~=15) and by the y value (12392->11637).
I hope my issue is understandable and that someone can help me.

(Sorry for my english)

Thanks

Raphaël


texts: Compare two points and conclude direction

texts:
Is there a simple method in geotools or JTS to conclude the direction between two geospatial points?

I would like to conclude if point A is (north, south, east, west, NW, NE, SW, SE) when compared with point B?

Or would I need to manually compute the angle between the points (possibly similar to http://docs.geotools.org/stable/userguide/library/referencing/calculator.html) and bin it into the aforementioned directions?


texts: Caching data locally using ArcGIS for Server?

texts:
Our organization currently works with our data over the network, which causes huge slow-downs.  We'd like to have everyone working off the data on their local drives to speed things up, and exploring ways to manage this.

I know ArcGIS Server should be able to cache feature and map tiles locally from a WFS/WMS but I am wondering if it will also allow an administrator to flag certain datasets to be cached to local machines so everyone's not always hammering the network each time ArcMap has to redraw the screen.

Is this possible?  If so, can an administrator make changes to the master networked data and have that change pushed out to all the local caches?

Basically I'm looking for some kind of source-control-like functionality which will allow people to work on local data, but be synchronized to a central dataset repository.  Selling my organization on moving into source control software is likely not going to happen, but they seem open to ArcGIS Server if that will solve the issue.

Any advice on how we can accomplish this?


texts: While reprojection a raster using 'Raster - Reprojection - Warp' in QGIS I can choose one of the following scanning techniques: Nah, Bilinear, Cubic, Cubic Spline, Lanczos. What does Nah mean and which one is recommended?


texts: Feature Class to Shapefile with definitions not values

texts:
I have a point feature class in an SQL Server Geodatabase that we use for day-to-day analysis. Occasionally we receive request for copies of our data from the public and send them shapefiles. If I export to a shapefile directly the coded values are what is written. I want to find a way to produce the shapefile with the textual description of the coded values in the coded value domain in the geodatabase.

  Shapefile:
  class = 00053

  Attribute Table:
  class = Gravel Street

It looks like ESRI's Production Mapping Extension has a tool that could do the trick but want to avoid purchasing an extension if possible.

How can I do this?


texts: ST_ClusterDBSCAN in PostGIS gives me 0 for the cluster ID

texts:
I am trying to cluster latitudes and longitudes from a dataset containing locations in Germany (n=52612, about 9000 without duplicates).
I am running the DBSCAN algorithm, but for none of the parameter values I have tried, I am getting any clusters- the cluster ID is 0 for all rows (I have also tried all coordinates vs. only unique ones).
I did a similar exercise with Python's DBSCAN in sklearn and got some clusters there for a broad range of parameter values.
What is wrong here?

eps = 20000 corresponds to 20km, is that correct?
Anyway, I don't get any results irrespective of the eps-value

--add geometry
ALTER TABLE table ADD COLUMN geom GEOMETRY;
UPDATE table SET geom = (ST_SetSRID(ST_MakePoint(lng,lat),4326)) ;
CREATE INDEX gix ON table USING GIST (geom);
--cluster
SELECT  *, ST_ClusterDBSCAN(geom::geometry, eps := 20000, minpoints := 100) over () AS cluster_id
FROM table
WHERE year = 2010
and country = 'Germany'


texts: How do I get the legend from a ArcGISDynamicMapServiceLayer using esri javascript api?

texts:
New Code: I updated the code a bit to more a recent recent version of the api, but ashima was correct in that I was not passing my layer as an array. Fixed that and it worked

 map = new Map("map", {
   basemap: "streets",
   center: [ -84.32, 46.5], // long, lat
    zoom: 13
});

//fetches myMap layer
var myMap = new FeatureLayer(myMapLocation, {
  mode: FeatureLayer.MODE_ONDEMAND,
  outFields:["*"]
});

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

map.addLayers([myMap]);

So I've been using the ESRI Javascript API, so far what I've been able to do is create a map, fetch a layer from another source and combine them, but I'm stuck on grabbing the legend from the other source and also adding that. This is what I currently have after my attempts (also tried switching the map : map to map : mymap)
The myMapLocation is just set to "/MapServer" and not "/mapServer/layers" or anything

 map = new Map("mapDiv", {
                    basemap: "streets"
                });
myMap = new esri.layers.ArcGISDynamicMapServiceLayer(myMapLocation);

map.on('layers-add-result', function () {
//this section breaks v
    var legend = new Legend({
        map: map,
        layerInfos: legendLayers
    }, "legendDiv");
// ^
    legend.startup();
});

legendLayers.push({ layer: myMap });

map.addLayers([myMap]);


texts: Is there tool to save currently visible WMS layer to file in QGIS, or some place in QGIS where you can get needed xml to fetch it using gdal ?


texts: Point-in-polygon using WGS84 geographic coordinate system

texts: When running a point-in-polygon operation, is it valid to run that with polygon features in straight up WGS84 and a point in lat/lng, or should both be projected into a non degree based coordinate system?


texts: OpenLayers plugin and QGIS 1.9: Error after install

texts:
After I installed the OpenLayers Plugin (1.0.0) on my QGIS 1.9 this error appears:

Konnte Erweiterung openlayers_plugin nicht laden wegen eines Fehler beim Aufruf seiner initGui()-Methode

Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 198, in startPlugin
    plugins[packageName].initGui()
  File "/home/map/.qgis//python/plugins/openlayers_plugin/openlayers_plugin.py", line 160, in initGui
    if not self.__setCoordRSGoogle():
  File "/home/map/.qgis//python/plugins/openlayers_plugin/openlayers_plugin.py", line 211, in __setCoordRSGoogle
    if not self.__coordRSGoogle.createFromEpsg(idEpsgRSGoogle):
AttributeError: 'QgsCoordinateReferenceSystem' object has no attribute 'createFromEpsg'

Python-Version:
2.7.3 (default, Sep 26 2012, 21:53:55)
[GCC 4.7.2]

QGIS-Version:
1.9.0-Master Master, exported

Python-Pfad: ['/home/map/.qgis//python/plugins/sextante', '/usr/share/qgis/python', '/home/map/.qgis//python', '/home/map/.qgis//python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-linux2', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PIL', '/usr/lib/python2.7/dist-packages/gst-0.10', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/pymodules/python2.7', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client', '/usr/lib/python2.7/dist-packages/wx-2.8-gtk2-unicode', '/home/map/.qgis/python/plugins/qgis2kml', '/home/map/.qgis/python/plugins/qgis2kml', '/usr/share/qgis/python', '/usr/share/qgis/python/plugins/fTools/tools']

Any ideas?


texts: Covisibility with QGIS/GRASS : r.los issue

texts:
I try to use r.los for a viewshed calculation.

I specify the dem (Aster DEM clipped on a 2km*2km region).

I specify x and y coordinates and run the command.

Nothing append. I waited for 20 minutes without a result. The loading cursor stay at 0. No error message, no stop, nothing.

I'm using r.los in QGIS/GRASS with an 1.8 version and under Windows 7.

Do you kwno why ?


texts: Adjusting raster size for usage in QGIS raster calculator?

texts:
I want to make a visibility analysis of planned wind turbines in my region.

I have a DEM raster and a polygon describing the possible area where the wind turbines might be set up.

People are complaining, that the view to a prominent mountain in the region will be disturbed.

My "simple" plan is, to raise the DEM-values covered by the polygon up to 200 meters and then perform a visibility analysis from the top of that mountain to be able to identify the regions, where the wind turbines might "block" the view.

Problem: I'm not able to successfully combine the vector and the DEM.

Steps taken so far:

Polygon was clipped from the DEM
Raster values from the clipped part were set to 200

What doesn't work:

When using the raster calculator, trying to add DEM + clip, the resulting image only covers the clip.

As i understand, the extend of the two rasters should be the same, and that NoData values will prevent operations (googling for 3 hours without success).

I fail in resizing my clip to the DEM and also do not get rid of the NoData values.(Really tried every tool under Raster to no avail)

Unfortunately, the rastercalc tutorials i found do not cover these "basic" things, so i'm completely lost at the moment.

I am using latest qgis master from osgeo4w.


texts: Can't find the ArcGIS License Manager folder

texts:
I am using ArcGIS 10.7 for desktop and my operating system is windows 10
I am trying to find out when my ArcGIS license expires, and this expiration date is supposed to be in the ArcGIS License manager folder. I can't find this folder when I am in the ArcGIS folder:
C:\Program Files (x86)\ArcGIS\Desktop10.7
There is just a License folder in the Desktop10.7 folder and it doesn't show the license information


texts: Assigning Sink ID to upstream features in ArcGIS geometric network?

texts:
Looking for a way to automatically assign an 8-character Sink ID to every upstream feature in a geometric network. There are over 500 sinks and about 5000 points and lines upstream from those sinks.  Would be great to click on any point or line and see what the Sink ID is, or quickly select records based on which sink they are upstream from.

A similar question is here:  How to assign attributes from hydrologic sinks to upstream features?


texts:
I have generated a virtual raster out of a large number of tifs but I am only interested in an area so I want to clip it. When I do it, it saves a raster file of several Gb. I would like to clip the virtual raster and obtain another virtual raster of the clipped area.

Is that possible?


texts:
Suddenly in my QGIS 2.x installations the official QGIS repository is unavailable from my computer. Also I can't load any WMS (from OpenLayers plugin etc.). It worked before and I don't know what I did that causes the problem.

No difference between OSGeo4W and stand alone installations.

I deleted QGIS in the registry and reinstalled QGIS but the problem persists.
I adjusted the firewall but the problem persists.

With QGIS 1.8 on the same machine there is no problem.

From other computers in the office there is no problem connecting to the repository and getting WMS tiles with QGIS 2.x (no proxy necessary, no firewall adjustments).

Any ideas?


texts: Create Voronoi/Thiessen polygons in QGIS around 3 points

texts:
Running QGIS 1.8.0
I have three bus compounds.  I would like to generate Voronoi/Thiessen polygons around those three points.  Following Voronoi/Thiessen adding an all school points will allow us to determine which school should be served by which bus compound.
Creating Voronoi:

It works very quickly and results in:

Obviously not what I expected.  Can anyone shed some insight?
Thanks


texts: Setup shadedrelief background in basemap using a cropped orthographic projection

texts:
I'm trying to plot a sector of the world in basemap using an orthographic projection and I would like to have a shadedrelief background.

If I don't use shadedrelief everything works fine:

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)

m = Basemap(projection='ortho',lon_0 = 0, lat_0 = 40,
        llcrnrx=-3000000, llcrnry=1000000, urcrnrx=3000000, urcrnry=7000000,
        resolution = 'l')

m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')

# draw parallels and meridians.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,420.,60.))
m.drawmapboundary(fill_color='aqua')

However, when I add the command m.shadedrelief() (after removing the m.drawmapboundary and the m.fillcontinents commands) I get the following graph:

Apparently the problem is that the background is not restricted to the bounding box I set in the basemap command. I tried also with other backgrounds (e.g. bluemarble) but they have the same problem.

Does anyone know how to fix this? I searched online but I didn't found anything unfortunately.


texts: PostGIS update statement returns wrong information

texts:
I have two tables in postgis with polygon data. Both are in wgs84. I want to update a column (column1) in the one table (table1) with information from a column (column2) in the other table (table2), where polygons in table1 intersects polygons in table2.

Here is the SQL:

UPDATE schema1.table1
SET column1 = column2
FROM schema2.table2
WHERE column2 IN (SELECT column2 FROM schema1.table1,schema2.table2 WHERE ST_Intersects(table1.geom,table2.geom));

The query runs, but it doesn't return the correct information. All the information in column1 is the same even though there are different entries in column2. I have also checked that the polygons in table1 intersects different polygons in table2, so the results in column1 after running the query is incorrect. I have tried this with intersects, overlap, within and contains with the same result.

What am I missing?


texts: Calculating distance to points on a coastline from all cells in extent

texts:
I'm on ArcMap 10.4.0.5524

Layers:

Layer 1: I have an outline (polygon) of an irregularly shaped island approximately 10km square in size.
Layer 2: I have a layer with various points snapped to the island polygon representing locations I want to calculate distance to, let's call them coastal towns.
Layer 3: I have the points layer (Layer 2) converted to raster. The extent of this layer has been changed to be the extent of a different layer approximately 200km square, which is all considered open sea. When I change NoData to be shown as a colour, this extent is displayed correctly. But when I zoom to layer, it zooms to where the coastal towns are instead of the extent of the whole layer.

The problem: I want to calculate the distance to the nearest coastal town from all grid cells in the 'sea'. I've tried to do this using the Euclidean Distance tool, which is probably the wrong one. What I get as output is a distance raster showing distances within the square expressed by the coastal towns furthest to the north, east, south and west, including areas on both 'sea' and 'land'. The maximum distance value returned is 2.56, which isn't correct and I want the values in kilometres.

How do I generate a distance raster using only routes available by sea (not crossing land, so I can't go with this and mask the result!), with results in kilometres and for the full 200km extent around the island instead of just the extent of the coastal town points?


texts: Transforming Data using 2D Costum Transformation in Arcmap 10.2 for Desktop

texts:
I have ArcMap layers of roads, buildings etc. in one reference system (Israel_TM).
I need to transfrom all this data (while keeping the topology) to pixel-space "reference system".
I have already calculated the transformation parameters (using Matlab) of a Conformal Transformation (2D) between the two reference systems: Israel_TM to Pixels. It it a 4 parameters transformation I have: 1xRotation, 1xScale, 2xTranslations.
How can I create a user-defined 2D transformation and use it to transform all my data to Pixel-space ?
All I could find is Helmert 7-Par. transformation or other 3D transformation methods under in Layer properties-->Coordinate System-->Transformation.
I tried to use the Toolbox for this but couldn't find anything relevant.


texts: Splitting line layer at polygon's boundary of polygon layer using PostGIS

texts: I have a line layer and a polygon layer in PostgreSQL 9.2, PostGIS 2.0. How can I split every line each time it crosses a polygon's boundary?


texts:
I have a census tract shapefile with the population of a minority group and total population.

I am trying to analyze this group by high and low risk flood zones.

The problem is that not all of the census tracts are completely covered by the flood shapefiles so the estimate wouldn't be accurate.

Is there a way to get an accurate estimate of the people in the flood zones?

Most of these tracts are urbanized.


texts: Should ST_OrientedEnvelope results differ depending on order of transform?

texts:
I have a PostgreSQL db with a table with geometry multipolygons stored in srid 4326 (lat/lngs).  I need to calculate the area of the oriented envelope around each polygon, and I am using st_orientedenvelope to do so.  I want area in reasonable units, so I also need to use st_transform before calling st_area

Why do I get different results depending on the order of transforming and generating the envelope?

select
    st_area(st_orientedenvelope(st_transform(geom, 26910))) a,
    st_area(st_transform(st_orientedenvelope(geom), 26910)) reverse
from
    mytable

Whether the envelope is generated and then transformed or the transformation happens and then the envelope is generated shouldn't (to my mind) matter.  Calculating a / reverse should always result in a number very close to 1.  However, in my data set of ~2000 records, more than 10% of the area calculations (a compared to reverse) differ from each other by more than 5% - some by as much as 25%.

In the image, the black line is envelope generated after transform (a), and the orange line is transform after creating envelope (reverse)

In the above example, I use 26910, which is SRID for UTM 10N because the data is in Seattle area.  However, I get very similar results using other SRIDs, including 2926 (Washington State Plane) and 2163 (US National Atlas Equal Area).


texts: More efficient way to get the start and the end points of a line (arcpy)?

texts:
I am collecting start and end points for a large polyline feature class in a file geodatabase. I am using arcpy, ArcGIS 10.3, Python 2.7.3 64-bit and a 10.3 compressed file geodatabase with just one feature class (and no other objects).

One accesses the geometry properties with the SHAPE@ token which gives access to the geometry of the line (namely firstPoint and lastPoint).

with arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"],query) as search_cur:
        allLinesList = [[row[0],(row[1].firstPoint.X,row[1].firstPoint.Y),
                         (row[1].lastPoint.X,row[1].lastPoint.Y)]
                            for row in search_cur]

Reading 1 million line features (note that my lines consist of only two vertices, start and end) takes ca 5 minutes. I have tried using for loop iterating a cursor instead of comprehension and iterating features in chunks of several hundred thousands with no significant difference in time.

Do you know any other faster way to read XY of the start and end points of the lines?


texts: Mobile "GIS"/geolocalization application for popularisation of science?

texts:
I am looking for an mobile multi-platform (or at least Android) application which will allow me to share my spatial data.

On day to day bases I am working with QGIS in archaeology and with Losuc GIS in the field. I have developed some interpretaions (remote sensing based) which I would like to share withe the public. I am looking for an application which would allow me to put a single link (QR code) on the information board that is placed on the archaeological site and which would refer to a mobile platform which would use geolocalization and which would contain my data (vector and if possible some raster). It has to be easy, simple and user-friendly. I would not like to share the data separately as I think it would be to complicated and not user-friendly enough. I am looking for an solution which would not require IT skills from me.

I've tryed using Google Earth online (earth.google.com) but it seems like the date is not apearing on mobiel devices (I guess its still in development phase).

Does such a system exists or should rather seek help within IT framework to develop specially for my purpose?


texts: Do Not Filter vs Ignore Invalid Geometries

texts:
In QGIS 3, What's the difference between "Do not filter (better performance)" and "Ignore invalid geometries"?

Both can be found at Settings > Options > Processing > General

Both are options to ignore invalid geometries during spatial analysis.


texts: Copying string data from one SpatialPolygonsDataFrame to another if they do not perfectly overlap

texts:
I have a SpatialPolygonsDataFrame of European NUTS 2006 regions (1462 regions) from here. I have imported this shapefile using the readOGR function in R. For each of the polygons, the NUTS 3 code is stored in the dataframe.
Before 2006, the regions were coded differently (e.g. Sweden's SE0A2 became SE232). Now I would like to add the 2003 NUTS codes to my 2006 SpatialPolygonsDataFrame. However, I cannot just use over(x = NUTS_2003, y = NUTS_2006), since the spatial objects do not perfectly overlap. The NUTS 2003 shapefile only has 1433 regions.
Does someone have any idea of how to tackle this efficiently?
This website explains how the NUTS codes have changed over the years, but manually entering all the 2003 codes does not seem to be the most efficient option.


texts: Currently, I have the NDVI values, and I transformed it from [-1, 1] to [0, 255]. Now my question is how do I visualize either of these value? My first thought is to make a color map similar to the image below. Are there any "official/mainstream" way besides assigning arbitrary colors to values?


texts: Seek assistance with: SQL Server 2008 Query that takes values from a column in another table (or from a flat file)

texts:
I am given a list in MS Excel of 554 statistical units (called Primary Sampling Units or PSUs) from another department and have to query all 554 PSUs from a SQL Server 2008 Spatial database table (PSUs_2012) that contains over 6000 PSUs.

I want to avoid having to type each of the 554 PSU numbers in the WHERE clause and is wondering if there is a way that I can write a query that takes values from a specified column (e.g. PSU_Number) in the Excel/CSV table and return all the rows with the same PSU numbers in the PSUs_2012 database table. I can also import the (excel/csv) data into the same database but would still have the same problem of how I should go about querying the main table (PSUs_2012) using all the values in a column of another table. Note that PSU numbers are unique IDs and are in both data sets with identical column names.

Please note that I am new to databases and programming but willing to learn, fast. Solutions using PostgreSQL and SpatiaLite are also welcome.


texts: Total loss of information while switching between layout view and data view? Arcgis for desktop 10.5

texts:
I managed to replicate my problem doing the following:

I have my MXD in layout mode, I open it and everything is fine as long as I stick to that view. When I change to data view and close the map, when I reopen it I have lost all information in the project, both in layout anda data view.

When I say that I lost the info, I mean that the table of contents has the data frames showing but the layers are missing, except in the second data frame which is an index map.

Sometimes I get an error messsage saying : "Unable to save. Check to make sure you access to the specified file and that there is enough space on the storage device to hold your file" although the problem still happens even if I dont save it and than reopen it.

I am using a custom map template. Could it be something related to it?

I am using the default geodatabase and it has probably hundreds of files there already, since its a big project and I did´t know better at that time.

Does someone have any idea of what is happening?

This my first time working with the program, can be a really basic user error as far as I know.


texts: How to optimise the query in PostgreSQL for function ST_Difference?

texts:
PostgreSQL

I executed below query for 4 million of geometry polygon in PostgreSQL. It took unlimited time which I need to stop the execution after 2 hour.

 SELECT
      c.id,
     c.classification,
    ST_Difference(b.geom, c.geom)
    FROM
      public."city" c,
    public.b
    where
    c.classificationin ('Class I','Sub Metro City','Class II','Metro City','Class III','Mega City');

QGIS

The same layer I used in QGIS using Difference in Vector >> Geoprocessing Tools and it took hardly 3 minutes to clip the region.

Please help me out

  how can i optimize the query in PostgreSQL ?


texts: Showing a point shapefile with raster layers in geoserver

texts:
I'm new to geoserver, and I am fiddling around with the geoserver suite for Windows. I have two aerial images georeferenced and imported as GeoTIFF which worked out fine. I am also trying to import a point shapefile with two points in it, but it wont show. The CRS of all three layers are the same, and if I just preview the shapefile, it looks good.

You can see how my setup looks here:


texts:
I would like to create an animated map/video in QGIS illustrating the route of a voyage (goig/panning from one stop to another). Something like http://www.animaps.com.

I don't think it is possible with Time Manager or MMQGIS? Maybe I'm wrong. With OpenLayers it should be feasible but it seems to require a lot of coding in Javascript and you still have to create the video afterwards from the screen animation.


texts: Progmatically remove unwanted fields from shapefiles without using disk

texts:
I want to remove specific fields that I don't need from multiple shapefiles stored directly in memory (the raw file data stored in a variable).

Unfortunately I haven't found any way to do this without having to write those files on disk, loading them, removing the unwanted fields, then reloading the files to memory and deleting the copy on the disk. And due to the shapefiles sizes, it takes a lot of time that could be shortened.

I only know PyQGIS so far, but I'm not bounded to it.

I actually do the following (using disk):

to_keep = ['Layers', 'I want', 'To keep']
vl = QgsVectorLayer(path, "tmp", "ogr")
attrib = vlIn.attributeList()

for elem in to_keep:
     i = vl.dataProvider().fieldNameIndex(elem)
     del attrib[i]
vl.dataProvider().deleteAttributes(attrib)
vl.updateFields()
QgsVectorFileWriter.writeAsVectorFormat(vl, path, "utf-8", crs, "ESRI Shapefile")

Note : Shapefiles are downloaded from internet and are zipped.


texts: Generate a map image file from tiles and features

texts:
I am looking for a way to generate images (e.g. png) from a tile source with features added on top of the map, without user interaction or requiring a browser.

Currently, I am able to create a leaflet map and show features on that map. Now I want to create an image of that map. But I want to create hundreds of different images, so I don't want to use a browser but execute it on a server.

Is there a way to do this? BTW: I don't depend on leaflet, I just need to be able to draw a feature on a map (e.g. POIs or lines).


texts: How to import ESRI geodatabase into QGis? (GIS beginner's question...)

texts:
Possible Duplicate:
  How to get “.GDB” (Esri File Geodatabase) support in Quantum GIS (osgeo4w - qgis)

First, I would like to apologize as I know that this topic has already got a lot of questions/answers in the last months. But despite those answers, I am stuck with my problem. The main reason might be I am not a GIS specialist and I don't have any computer programming background.

I've received an ESRI geodatabase and I would like to open it with QGis. I tried to open it via the "Add layer" checking the "Directory" box but I always get the answer "It's not a valid or recognized datasource".

Reading the previous comments I assume I have to set up something, QGis is not able to read ESRI geodatabase with its native configuration.

Does anyone could give me any guidance?

I have QGis 1.8.0 Lisboa that I installed through OSGeo4W using the Express Deskop install option. I have Windows 7.

Thank you

Best regards,

Bap


texts: What is the easiest way to install libgdal 3.3.0 and header files (gdal-devel) on Ubuntu 18.04?

texts:
I would like to pip install GDAL 3.3.0 Python package, but it requires libgdal (3.3.0 or greater) and header files (gdal-devel). APT package manager repositories only have gdal-bin 2.4.2 and libgdal-dev 2.4.2.
Is there any other way to install libgdal 3.3.0 than from source?
If installing from source is the only option, is there any step-by-step guide for that? I once tried it by installing sqlite3 which was required by PROJ 6 and then PROJ 6 which was required by GDAL and still got configure: error: PROJ 6 symbols not found error. I guess I didn't just understand something.


texts: I am a beginner to QGIS. My problem is I tried to create sub-catchment within a basin by using provided point datasets (river gauge location). In total, I have 12 measuring points (i.e 12 pairs of x and y coordinates) within my basin. I want to delineate all catchments within basin on the basis of provided coordinates ( not on the basis of minimum or maximum cell that r.watershed has.) I am looking for the function in 'hydrology' extension tool of ArcMap.


texts:
I want to perform Slope using PostGIS. I used this query but I got a strange result unit:

create table my_slope as
SELECT ST_Slope(rast,1,'32BF','PERCENT',1, false) as slope
FROM my_raster

I wonder know why I got this strange result.


texts: Issue while running grass from subprocess

texts:
I have a map based web application using geodjango. I needed to display the crime hotspots within a selected polygon area. I found that v.kernel is very good for generating heatmaps/hotspots. So in my application I ran this command through python subprocess and created the raster map showing point density. I displayed this raster map over the polygon area in grass. But when I try to run the command through python using d.out.file command, error is raised.

command5 = "d.mon start=x1"
command6 = "d.mon select=x1"
command7 = "d.mon unlock=x1"
command8 = "r.colors byr"
command9 = "d.rast crimes_"+str(pname)
command10 = "d.vect "+str(pname)+" type=line,boundary"
command11 = "d.out.file output="+str(pname)+" format=png size=800,600"
command12 = "d.mon stop=x1"

saveCommands = [command4,command5,command6,command7]

for sc in saveCommands:
    s = subprocess.Popen(sc,shell=True,env=environ)
    s.wait()

Error is :

ERROR: Graph_Set: can't open Display
ERROR: Unable to get GIS_LOCK enviroment variable value
Problem selecting x1. Will try once more
ERROR: Unable to get GIS_LOCK enviroment variable value
d.mon select=x1
Exception AttributeError: AttributeError("'_DummyThread' object has no attribute '_Thread__block'",) in <module 'threading' from '/usr/lib/python2.7/threading.pyc'> ignored
ERROR: Unable to get GIS_LOCK enviroment variable value

All I need is to export the grass vector and raster layer to an image.


texts: When building fieldmappings, how do I handle the SHAPE field?

texts:
I have a collection of hundreds of shape files and a python dictionary mapping the input field to the output field. I'm building my fieldmappings dynamically using this dictionary since it already exists and is in a configuration file, editable by the user if fields change.

How should I handle the SHAPE field? I want the geometry copied over during any merge operations. The mapping in my config file, "SHAPE@": "SHAPE@" (input field : output field) throws an error for the shape field but all other mappings work.

Update 3: Code included...

field_names = GetFieldNamesInFeatureClass(file)

fieldmappings = arcpy.FieldMappings()

for input_name in field_names:

    output_name = fmFromCfg[input_name]

    fieldmap = arcpy.FieldMap()
    fieldmap.addInputField(file, input_name)

    fld = fieldmap.outputField
    fld.type = "DOUBLE"
    fld.name = output_name
    fld.aliasname = "Output Field Alias"
    fld.isNullable = True

    fieldmap.outputField = fld
    fieldmappings.addFieldMap(fieldmap)

The configuration file is JSON. It is read into fmFromCfg. Here is a simplified version with relevant section:

{
  "FieldMapping": {
    "FID": "OWNER",
    "SHAPE": "SHAPE"
  }
}

Update 2: Changed the field mapping to "SHAPE": "SHAPE" but I'm still getting an error, RuntimeError: FieldMap: Error in adding input field to field map.

Update 1: When I change the mapping to "SHAPE": "SHAPE@", a match is made but I get RuntimeError: FieldMap: Error in adding input field to field map because the @ at the end of the output field.


texts: Add raster layer or a WMS using geotools - problem to run the jar file from cmd

texts:
I am using geotools in order to develop a custom GIS.

I tried to add a raster layer and a WMS based on the tutorials given by geotools. When i am running the project using the main class in Netbeans or Eclipse everything is working properly. However, when i tried to run the jar file from the cmd, the above functionalities crashes. Below you can see a screenshot which shows the error regarding the WMS.

I forgot to send that the project which i am working on is a Maven project.

Can you help me to understand where is the problem?

<dependencies>
    <dependency>
        <groupId>junit</groupId>
        <artifactId>junit</artifactId>
        <version>4.11</version>
        <scope>test</scope>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-main</artifactId>
        <version>13.5</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-shapefile</artifactId>
        <version>13.5</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-swing</artifactId>
        <version>13.5</version>
    </dependency>
    <dependency>
    <groupId>org.geotools.xsd</groupId>
        <artifactId>gt-xsd-core</artifactId>
        <version>9.3</version>
    </dependency>
    <dependency>
        <groupId>org.geotools.xsd</groupId>
        <artifactId>gt-xsd-kml</artifactId>
        <version>9.3</version>
    </dependency>
    <dependency>
          <groupId>org.geotools.jdbc</groupId>
          <artifactId>gt-jdbc-postgis</artifactId>
          <version>13.5</version>
    </dependency>
    <dependency>
      <groupId>org.geotools</groupId>
      <artifactId>gt-geometry</artifactId>
      <version>${geotools.version}</version>
    </dependency>
    <dependency>
    <groupId>org.postgresql</groupId>
        <artifactId>postgresql</artifactId>
        <version>9.3-1103-jdbc41</version>
    </dependency>
            <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-epsg-hsql</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-geotiff</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-image</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-wms</artifactId>
        <version>${geotools.version}</version>
    </dependency>
</dependencies>


texts:
I have a large file of a road network. I need to calculate the value of direction (in degrees or radians) of the line segments, and save this value to a new field. I need to use MapInfo.

I have found this excellent example using qgis, but need to undertake the work using MapInfo. Does anyone have a solution. I tried the basic idea of using the function listed in the solution and copy and pasting that in, but it does not work.

How to add Direction and Distance to attribute table?

Hope this all makes sense. To conclude, I need to calculate direction for each line segment and add them to a new field in MapInfo.


texts:
I'm unable to print a AT-Basemap-Layer with the GeoServer Print Plugin

I could not find any working example with WMTS-Sources.


texts: There are different types of spatial diffusion such as contagious model, migration model, etc.  What are the good books to learn about the mathematics behind these spatial diffusion models?


texts:
I have a few dozen feature classes, one of which contains the field I am looking for. They are stored in several different file geodatabases.

Is there a quick way of searching each feature class to find the one that has my field? Or do I have to check each one's attribute table in Catalog? Some of the files have the same names but are stored in separate locations.

I would like to be able to do this inside a Calculate Value tool in ModelBuilder, so that each time I needed to, I could run the model on a single file gdb and it would search through all the feature classes within that gdb and look for the field.

So far, I have a feature class iterator which iterates over all the feature classes in the file gdb recursively. Inside the Calculate Value I have the following code:

findField(r"%Feature%", "%Search Field%")

def findField(fc, fi):
  lst = arcpy.ListFields(fc)
  for f in lst:
    if f.name == fi:
      return fc
    else:
      return "not found"

The returned values are collected in Output Values as strings. The output for all the feature classes has been "not found", even though the field exists in at least 4 of them.


texts: Intersect of all polygons from the same layer in QGIS

texts:
I have various polygons within one layer in QGIS.

As a result I need the area where all polygons overlap as shown in the picture.
Is there a way to do?


texts: Does anybody have a link to a video or tutorial on how to georeference an image file with a crossection in ArcScene?

texts: I am looking to put together a 3D site conceptual model. I have subsurface data as well as surface data and would like to display it in ArcScene. How do I georeference a crossection so that it displays in 3D?


texts: Bounding box along parallels and meridians

texts:
I have a PostGIS table with a column location geography(Point,4326). To select data from a specific region I use location && ST_MakeEnvelope(65,-77,120,-75).
This may work well for most applications, but my data is in the Antarctic region and there it is not exactly doing what I want.
What the intersection does is creating a rectangle with corners according the boundaries as you can see in the image attached in red. What I want is a selection of all points between latitude -75° and -77° and longitude 65° and 120° (the black polygon).
Is there any possibility in PostGIS to do such a query where the bounding box isn't a rectangle but follows the parallels and meridians?

EDIT

Having a closer look at the documentation of "ST_MakeEnvelope" and trying it myself seems like ST_AsText(ST_MakeEnvelope(65, -77, 120, -75, 4326)) simply translates to "POLYGON((65 -77,65 -75,120 -75,120 -77,65 -77))". The corners of the polygon are logically connected by the shortest line on the sphere/spheroid which is a great circle/geodesic. This is exactly the behaviour I would expect for any test with a polygon, but it's not what I expect from a test for lonmin, latmin, lonmax, latmax.
It seems like ST_MakeEnvelope is not the correct function for such a query, especially close to the geographic poles.
Is there any other possibility without converting the column to geometry (which probably prevents the usage of the index on this column and lasts way too long)?


texts:
Is it possible to draw one way road direction arrows using CSS of GeoServer?

I want to get the same result as shown below (see img-1).

(UPDATE)
So far, I've added arrow sign as label in separate layer and this is my result (see img-2). Arrows appear not everywhere, they are missing on "Guramishvili st" for instance.


texts: Calculating mean for each buffer using QGIS?

texts:
I am attempting to calculate mean NDVI at different time points with 100m and 500m radii of certain co-ordinates. Thus far I have used 'Raster Calculator' to calculate the NDVIs and have used 'Vector -> Geoprocessing Tools -> Fixed distance buffer' to calculate the buffers. My screen looks like:

How do I calculate the mean NDVI for each layer in each buffer?

I am hoping to end up with a 9 column csv table I can then import into R with the column headings:

  id; lat; long; ndvi_100m_summer1990; ndvi_500m_summer1990;
  ndvi_100m_spring1990; ndvi_500m_spring1990; ndvi_100m_summer1991;
  ndvi_500m_summer1991


texts: Splitting service area polygons in ArcGIS Pro

texts:
I am doing a project in which I am looking at the areas accessible within varying distances of open subway entrances. When I created the Service Area layer for open subway entrances, I noticed some peculiar things. The polygons were designed to be split, or not overlapping.

I set cutoffs for 1324 feet, 2728 feet, and 4092 feet from the facilities, which are open subway entrances. There should be polygons for 0-1364 feet, 1364-2728 feet, and 2728-4092 feet. For the most part, this is what has happened, but in some areas, there are weird errors.
For instance, at Union Turnpike station, the polygon for the areas within 0-1364 feet of the 80th Road subway entrance, has 0 area, and I am instead in the 1364-2728 polygon.

Another error is that though I set the service polygons to be split, they are overlapping, which is incorrect.

There are also issues like this:
For the same station entrance, you have multiple overlapping polygons. For instance this entrance has polygons covering 0-1364 feet, as it supposed to, and for some reason, one from 0-1273.15. Even weirder, these two areas are very different areas.

There also several disconnected service areas that don’t make any sense. On the whole, these errors might not make much of a difference, but you never know.

This one is also very peculiar. It should not be possible for a service area layer to include more than one facility, but this one covers two facilities, which is incorrect.

Does anyone have any idea of what is going wrong/what to do to fix the issue?
Edit: I tried using the detailed service area polygons and used lines + polygons, and it looks a lot better.

I have been adjusting the polygon trim distance to see how that changes things, and it does make things more accurate. However, there is an issue.
Areas in the center of blocks, which are not directly accessible by the street network show up beyond the threshold they are supposed to. As such, the centroids of the census blocks will be in the wrong polygon. Is there any way to merge segments like these in the middle of blocks to the polygons that surround them without doing it manually, one by one?
Also, if you look at the two images below, you will notice something weird. The two purple circles are subway entrances-one at the northeast corner of Nostrand Avenue and President Street, and one at the northwest corner of that intersection. For whatever reason, areas on the east side of Nostrand Avenue (the north-south street in between the two entrances), are shown as part of the polygon for the northwest corner, which does not make any sense, as it would take less time for people on the east side of the street to get to the entrance at the northeast corner. Does anyone have any idea of what is going wrong/what to do to fix the issue?


texts: Algorithm to implement for batch processing in ArcMAP modelbuilder

texts: I have a directory with about 10 polygon shapefiles. I need ArcMAP to take all these shapefiles and union them all into one layer using the model builder. I know I need to use some sort of iterator, however the basic one does not successfully union them all, but rather only the last one. What am I doing wrong here? Here is the model:


texts: Is there free pipeline shapefile data for Asia?

texts: I need to find oil & gas pipeline data in the world especially about Asia but I want to get them freely.


texts: ValueError using MapProxy decorate image middleware with WMS

texts:
I am successfully using Mapproxys decorate image middleware to apply a watermark to WMTS tiles. I am using the following code:

def annotate_img(image, service, layers, environ, query_extent, **kw):

    # Required opacity of watermark
    opacity = float(0.5)

    # Get path to watermark png
    wm_path = Image.open(os.path.join(os.path.dirname(__file__), 'watermark_256.png'))

    # Grab tile and Convert to RGBA
    img = image.as_image()
    if img.mode != 'RGBA':
        img = img.convert('RGBA')

    # Create new RGBA image for watermark
    wm = Image.new('RGBA', img.size, (0,0,0,0))
    wm.paste(wm_path)

    # Adjust opacity of watermark
    alpha = wm.split()[3]
    alpha_c = ImageEnhance.Brightness(alpha).enhance(opacity)
    wm.putalpha(alpha_c)

    # Create new composite image using the wm as a mask
    compo = Image.composite(wm, img, wm)
    final_img = Image.new('RGBA', img.size, (0,0,0,0))
    final_img.paste(compo)

    return ImageSource(final_img, image.image_opts)

As I said this is working fine for WMTS, however, when I try to view a WMS layer I get the following error in my log:

self.im.paste(im, box)
ValueError: images do not match


texts: Warping Image having Lat Lon Geolocation array

texts:
I'm trying to implement a script to project an image having the two matrices of Lat and Lon points. I've followed the steps proposed in:
Unable to warp HDF5 files
and everything is fine. I create the three .vrt files and than by using:

gdalwarp -geoloc -t_srs EPSG:4326 img.vrt output.tif

the final image is correctly projected.
Also, just for visualization purposes the following code does is job:

plt.pcolormesh(lon,lat,img, cmap='rainbow')
plt.show()

Now, I want to create a script without having to create the .vrt files every time. So I started following Writing numpy array to raster file
but I'm not able to obtain the desired projection (Actually, the output image is similar to the input one, just horizontally translate of some pixels). Here is the key part of code:

ny = img.shape[0]
nx = img.shape[1]

xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()]
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymin, 0, -yres)

# create the 1-band raster file
dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF.tif', ny, nx, 1, gdal.GDT_Byte)

dst_ds.SetGeoTransform(geotransform)    # specify coords
srs = osr.SpatialReference()            # establish encoding
srs.ImportFromEPSG(4326)                # WGS84 lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file

dst_ds.GetRasterBand(1).WriteArray(img)   # write r-band to the raster
dst_ds.FlushCache()                     # write to disk

dst_ds = None

I think I'm missing something when I create the geotransform matrix because in this way I introduce only one corner of the image without giving any info about the transformation.

How do I solve the problem?


texts: How to do regression analysis in SAGA for topographic correction

texts:
Currently I am trying to do the topographic correction using C-Correction. I would like to know if there is a way to do the regression analysis between L (uncorrected reflectance of each band from satellite image) and cos i (the incidence angle), to get the C value.

I tried to do the regression analysis using SAGA (the uncorrected reflectance of each band as dependent variable). but I am not sure how to calculate cos i as the independent variable.
Do you have any advice how to get the "cos i" value?


texts: ArcObjects ITextSymbol.GetTextSize in Python revisited: usage with StartDrawing and FinishDrawing methods

texts:
I am using ArcGIS 10.0/ArcObjects and Python 2.6.  I am trying to get the usage of the ITextSymbol.GetTextSize method correct.  The syntax, as established in an earlier question, is:

dblx, dbly = pTextSymbol.GetTextSize(pDisplay.hDC, pDisplay.DisplayTransformation, txt)

and the outparameters are returned just fine.  My problem is in figuring out how to use this method so that it does not screw up ArcMap's screen cache.  ESRI documentation indicates that it should be surrounded by calls to IDisplay.StartDrawing and IDisplay.FinishDrawing so that the input hDC will be correct.  (Indeed, calling GetTextSize without the start and finish drawing results in an invalid hDC.)  The syntax that I am using for StartDrawing is as follows:

pDisplay.StartDrawing(0, esriDisplay.esriNoScreenCache)

I have tried both 0 (as discussed in documentation) and pDipsplay.hDC for the first parameter of StartDrawing, but the ArcMap screen still gets stuck in one place.  (Panning just causes a redraw from cache right back to the same place.)  Either input for the first parameter results in valid outparams but a stuck screen.  The problem could also be with the esirNoScreenCache constant, I suppose, but I don't see how.  I am not trying to draw anything to the screen at this point.  I'm just operating on a TextSymbol, not a TextElement or the like.  pDisplay and txt are both valid, in-process objects, so the problem, I think, is with the hDC or the screen cache.  I'm at a loss to figure out which.  Any insights are greatly welcomed.


texts:
I'm pretty new to Ubuntu and cannot figure out how to resolve this issue. I'm running the latest Geoserver WAR (2.13.0) on Tomcat7 on Ubuntu 16.04.4. I've added the GDAL library in order to support ECW files on Geoserver. When i try to create a new layer in Geoserver with ECW, the server crashes. The tomcat log /var/log/tomcat7/catalina.out shows the following error:

terminate called after throwing an instance of 'std::length_error'
what():  basic_string::_S_create

When i run:

/usr/lib/javainfo/gdalinfo.sh --formats | grep ECW

as expected it shows the ECW support:

ECW (rw): ERDAS Compressed Wavelets (SDK 3.x)
JP2ECW (rw+v): ERDAS JPEG2000 (SDK 3.x)

When i test it with a sample ECW file:

sh /usr/lib/javainfo/gdalinfo.sh /home/ubuntu/TILES_CM_062.ecw

I get the error below:

terminate called after throwing an instance of 'std::length_error'
what():  basic_string::_S_create Aborted (core dumped)

I can successfully create a GeoTIFF layer and when i check a sample GeoTiff with /usr/lib/javainfo/gdalinfo.sh there is no issue.

When i run:

/usr/lib/javainfo/gdalinfo.sh --version

it shows:

GDAL 1.9.2, released 2012/10/08

Here's the background to what i've tried so far. I followed these instructions at https://geoserver.geo-solutions.it/edu/en/install_run/imageio_ext_install.html. Specifically I have done the following:

Downloaded https://demo.geo-solutions.it/share/github/imageio-ext/releases/1.1.X/1.1.16/native/gdal/gdal-data.zip
Extracted into /usr/local/share/gdal/data
Set the GDAL_DATA environment variable (GDAL_DATA=/usr/local/share/gdal)
Also set GDAL_SKIP=JP2ECW as per https://demo.geo-solutions.it/share/github/imageio-ext/releases/1.1.X/1.1.16/native/gdal/linux/notes_on_ECW_and_JP2ECW.txt
Downloaded https://demo.geo-solutions.it/share/github/imageio-ext/releases/1.1.X/1.1.16/native/gdal/linux/gdal192-Ubuntu12-gcc4.6.3-x86_64.tar.gz
Extracted into /usr/lib

I also went through http://docs.geoserver.org/stable/en/user/data/raster/gdal.html.

I also tried installing the latest GDAL:

 sudo apt-get install gdal-bin

Running this again shows:

gdal-bin is already the newest version (2.1.3+dfsg-1~xenial2).

I think the issue must be to do with the GDAL library install?


texts: Show distinct features count as label in ArcMap10.2

texts:
How to show features count and not the value as a label in ArcMap10.2?

I have a polygon layer which contains 100 or so features. All the features are grouped by a group_name field containing values like "Forest", "Park", "Factory", etc.

I want to show the number of records present under a group. Suppose two features contains "Park", four features contain "Forest" in the group_name, so the label should be 2 on those features which contain "Park" and 4 on those features which contains "Forest".

EDIT:
The approach should be dynamic, so that, if some feature added/deleted then layer must reflect the difference. Like, if someone add an other feature and added a group_name= "Park" (in Attribute Table). Then the label should be automatically show 3 on map. I'm not concerned with TOC as such because I want to set the label expression and then I will publish it. The application users can add or delete the feature and they want to see the number of records as feature's label.


texts: Object height determination through Split Stereo image visualization QGIS

texts:
Is there any way to display (Aero-triangulated) left and right images (Stereo photogrammetry) in 2 mapview/viewport(different viewport/Map View/Multiple Canvas) in QGIS. So that we can use Sensor parameter (Azimuth and Elevation angle of different satellite) to calculate Object height By syncing of Left and right image.

Clicking on Point A (Left Image) and same point on Right image will calculate Azimuth from top to bottom of wall for both image is height of Object (Please correct if I am wrong)
Multiplication factor also can be applied


texts: Haversine distance versus Euclidean on an eqc "equi-distance" projection

texts:
I've got a network covering a large area, but the individual links are fairly small (<1km).  For calculating edge lengths I'm trying to decide whether it would be better to use Haversine distance on the decimal degrees or Euclidean distance on all the geometries converted into a CRS designed for distance measurements.
Option 1: Haversine Distance on the (lon, lat) of endpoints in 'epsg:4326' (python code for reference):
####==== Calculate the great circle distance in meters for two lat/lon points
def haversineDist(lon1, lat1, lon2, lat2):
    # convert decimal degrees to radians
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    theAngle = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    return 6367000 * 2 * math.asin(math.sqrt(theAngle))   ## distance in meters

Option 2: Pick a center point to set the CRS, convert the geometries to that CRS, and calculate edge lengths using Euclidean distance. For example, use the CRS +proj=eqc +lat_0=35.6812 +lon_0=139.7671 +units=m for the area around Tokyo.
I think Option 2 is more accurate for areas fairly close to the reference point, but what if I am also measuring lengths of edges several degrees away, such as around 43.52, 141.62?  If I want to keep the measurements of length consistent and easily reproducible, then Option 1 seems better.
I am still fairly new to these considerations, so there may be even better options that I am not aware of.


texts: ModelBuilder - Iterate through tabels to join and export and keep original field names

texts:
I am trying to iterate through an excel with over 200 sheets and join each to a ZIP code layer. I then need to symbolize each of those the three different ways and save that to a layer file. In order to do that, I need the model to keep the original field names so the symbology can be done for all. Here is what my model looks like

I have tried multiple routes to keep the field names. Utilizing the field info variable from Make Feature Layer, I can have the fields from the layer I am joining to (ZIP 2014) keep their original names, but the fields from the table still use the "Sheet__fieldname" convention and the fields from the table are what I have to base the symbolization on. I have unchecked maintain fully qualified field names for the Model properties as well as for the individual tools, and that doesn't work.

My only workaround has been to export everything in the model as is, then create another model and use feature class to feature class to change the field names. I truly believe there has to be some way to get all this done in one swoop.

I am using ArcMap 10.2.2