texts:
I'm trying to remove inner polygons with postgis.
I've created a routine, that resumes in:

From a polygon table, with a unique id:

Create a linestring table from polygon
Import to postgis
Get the length of each line

Now I have a table with, for example, 5 linestrings for each id, each linestring was a polygon, and I want to get ONLY the longest line from each ID.

For this, I have a select, but I don't know how to group the geometry, to get ONLY the longest line

select  id, geom as length from esp_lines group by id

I have to group it by another field, but I don't know which one.

The other solution I've found is:

select  id, max(geom) as length from esp_lines group by id

but max(geom) doesn't give me the expected result (the longest linestring)

Does anybody know if is there any function to get the longest linestring of a group?


texts: Service Area: Multiple network locations per facility

texts:
I am generating service areas using ArcMap's Ready-To-Use Services. I'm Using ArcMap 10.3. I am doing walking distance to train stations using 500m 1000m and 1500m.

I manage to create the service areas, but I find the locations (train stations in this case) snap to only one point on the network. This results in the service area being created only on one side of the train line (the side where the road is closer) despite being able to access from both sides in reality (see example below).

Is there a way to make the location snap to more than one road on the network without adding points manually (I have many train and bus stations that I am doing this with)?


texts: How to recover from deleting the symbolid of a arcgis annotation feature class?

texts:
I have arcgis 10.1 (arcinfo license) and I am using a versionned database (arcsde) with SQL Server.  I have accidentally deleted the symbolid of one of my annotation feature classes.  All annotations using that symbolid now have a red hallow around them and appear empty (although the attribute table is OK and has all the data).  I am trying to recover from this problem.  The following solutions do not work:

use field calculator to assign a new symbolid to the annotation class.  If I do this, all annotations re-appear but at zero angle (data in the "angle" column is lost and replaced by zero).
select the problem annotations from by backup database (in which these annotations are ok) and append them to the actual database.  For a reason that I do not understand new columns are created in the actual database (the annotation class had several custom columns added to it at the moment of creation and they seem to be doubled in the append process).  This doesn't look like a promising solution.

My fallback plan is to use the backup database and replace my annotation layer with the old one, but that means loosing 1 days work which I want to avoid if this is possible.


texts: How to read data, https using 7001 port in a NAS drive with QGIS?

texts: I have some geograhical data in a Network Attached Storage. How can I read that using QGIS through https protocol by 7001 port ?


texts:
Having trouble getting my code to work. Keep getting IOError that my parks.shp doesn't exist....not sure I understand what this is saying.

import arcpy

mxd = "D:\...\TravisCountyAustinTxx.mxd"

mapdoc = arcpy.mapping.MapDocument(mxd)

dataset = "D:\...\TravisCountyAustinTxx.mxd\parks.shp"

spatialref = arcpy.Describe(dataset).spatialReference

extent = arcpy.Describe(dataset).extent

for df in arcpy.mapping.ListDataFrames(mapdoc):

      df.spatialReference = spatialref
      df.panToExtent()
      df.scale = 15000
      print "\nData frame " + df.name + " contains the following layers:"
      lyrlist = arcpy.mapping.ListLayers(mapdoc, "", "", df)
      for lyr in lyrlist:
          print lyr.names
mapdoc.save()
del mapdoc

Full error message:

File "D:\Python27\ArcGIS10.2\Lib\site-packages\Pythonwin\pywin\framework\scriptutils.py", line 325, in RunScript
    exec codeObject in main.dict
  File "D:\CLASSES\Debugging_andErrorHandlingLab\Scripts\Mod3_Script2_Template (2).py", line 12, in
    spatialref = arcpy.Describe(dataset).spatialReference
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy__init__.py", line 1234, in Describe
    return gp.describe(value)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\geoprocessing_base.py", line 374, in describe
    self._gp.Describe(*gp_fixargs(args, True)))
IOError: "D:\CLASSES\Debugging_andErrorHandlingLab\Data\TravisCountyAustinTxx.mxd\parks.shp" does not exist


texts: Divide polygon layer by point layer in QGIS?

texts:
I have two layers - one polygon layer that is a buffer of a road layer, and one point layer that shows certain locations along that road.

What I want to do is to divide/split the polygon layer at the exact locations of the point layer. I also need the divided polygon layer to include attributes from the point layer.

Any ideas on how to do this in QGIS?

EDIT: I need the buffer to be 5-15 meters from the road.

Thanks!


texts: 'Evaluation error' when using Aggregate tool in QGIS

texts:
In QGIS 3.16. I am using the Aggregate function to reduce the number of features in a layer. I am grouping by one field and I want to concatenate the others.
My actual data table is fairly large and the operation seems to be working for some of my fields, but for some fields, I get an error message related to the concatenation, and for the life of me I can't figure out why.
Here is a screenshot of my Aggregate setup for a truncated version of my table, with a few of the features that are having problems. I am grouping by the field "Complexe"; note that all output field types are set to String (Text) with a length of 100, so that is not the issue:

This is the message that I get :

Evaluation error in expression
"concatenate("ID_CH","Complexe",TRUE,',')": Could not calculate
aggregate for ID_CH. Execution failed after 0.86 seconds.

Here are the files if you'd like to try and reproduce the problem: https://drive.google.com/drive/folders/1nOV1dxZWzNiIg9F77brQ36Rcvoze8QSb?usp=sharing


texts: Finding locations of highest values in raster using ArcGIS Desktop?

texts:
Using ArcGIS 10, I have a raster where I would like to find the pixel with the maximum value in the raster and return its location (center of pixel) in decimal degrees. I would like to iterate through this process returning the location of the second highest value of the raster, then third, and so on so that in the end I have a list of N locations that have the highest values in the raster in order.

I imagine that this might be most easily done using a Python script but I'm open to other ideas if there is a better way.


texts: Using PostGIS Intersection to carry fields from input layers into output layer

texts:
I'm having trouble carrying over the attributes from the input layers in my PostGIS intersection below:

CREATE TABLE intersection_output AS
(SELECT ST_Intersection(a.the_geom, b.the_geom) FROM nyct2010 as a, test_buffer as b WHERE ST_Intersects(a.the_geom, b.the_geom))

I specifically want to carry over the field geoid from nyct2010 and observation from test_buffer.

I've tried this, but doesn't seem to work:

ST_Intersection(a.the_geom a.geoid, b.the_geom.observation)


texts:
i have just created Tiles with Qtiles on Qgis, but how can put these folders (tile .png) on GeoServer or Mapserver?

I am trying to put the folder in  C:\Program Files (x86)\GeoServer 2.6.0\data_dir\gwc but i can't see it on GeoServer.


texts:
In some maps, the orientation font in the toolbar get messy( I work with Window 7 professional):

I tried to open a new project after i restart the computer- still it as not solved the problem.
I added those 2 pictures for the example, but it occurred also in all other toolbar. Even when i change the style into window (in option--> general) the problem still remain.


texts: Python error with OpenLayers plugin QGIS 2

texts:
I've been using QGIS for a while but since using QGIS 2 the OpenLayers plugin doesn't work for me on my home PC. However, I have installed the same version at my workplace with success. I have looked through other posts, and tried setting the CRS and allowing on the fly but it doesn't work for me.

Can anyone offer any other guidance?

I have also tried uninstall/ re-install.

Error code below:

  An error has occurred while executing Python code:

  Traceback (most recent call last):
    File "C:\Users\Lee/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 48, in addLayer
      self.plugin.addLayer(self)
    File "C:\Users\Lee/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 194, in addLayer
      layer = OpenlayersLayer(self.iface, self.__coordRSGoogle, self.olLayerTypeRegistry)
    File "C:\Users\Lee/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py", line 64, in __init
      self.page = OLWebPage()
    File "C:\Users\Lee/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py", line 38, in init
      proxy = getProxy()
    File "C:\Users\Lee/.qgis2/python/plugins\openlayers_plugin\tools_network.py", line 42, in getProxy
      proxy.setPort(settings.value("/proxyPort", 10, type=int))
  TypeError: unable to convert a QVariant of type 10 to a QMetaType of type 2

  Python version:
  2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

  QGIS version:
  2.0.1-Dufour Dufour, d94c044

  Python path: ['C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing', 'C:\Users\Lee\.qgis2\python\plugins\GeoCoding', 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python', 'C:\Users\Lee/.qgis2/python', 'C:\Users\Lee/.qgis2/python/plugins', 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins', 'C:\PROGRA~1\QGISDU~1\bin\python27.zip', 'C:\PROGRA~1\QGISDU~1\apps\Python27\DLLs', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\plat-win', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISDU~1\bin', 'C:\PROGRA~1\QGISDU~1\apps\Python27', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg', 'C:\Users\Lee\.qgis2\python\plugins\mmqgis/forms', 'C:\PROGRA~1\QGISDU~1\apps\qgis\python\plugins\fTools\tools']


texts: Descartes Labs Platform scenes.search within less cloudy periods

texts:
Sometimes searching for scenes over a large date range is necessary to obtain cloud-free images, but can be prohibitively slow for satellites with high revisit rates (e.g. Sentinel-2).
Is there a simple way to search for scenes with scenes.search to only return scenes within less cloudy time periods (e.g. summer) at the location?


texts: What is the "57.29578" that ArcGIS uses in Aspect calculation?

texts:
ArcGIS Desktop Help on How Aspect works explains the calculations used, but introduces the for me magic number of 57.29578:

aspect = 57.29578 * atan2 ([dz/dy], -[dz/dx])

Where does this number come from and what does it stand for?


texts: How to get Google Street View panoramas on the same street as the requested address?

texts:
The Google Maps API geocodes an address to a point. I can use that point to find the nearest panorama (street view) image. I can then determine the heading based on the two points and show the panorama pointed at the right house.

The problem is that in some cases such as corner lots, the closest panorama might be on a different street. For example, a house with an address of 654 NW 4th Street is on the corner of 4th and Greenwood.  The closest panorama is showing on Greenwood instead of 4th so I get a picture of the side of the house instead of the house.

Am I doing something wrong here?  Is there a way to get a list of panorama's so that I can find one that is on the correct street?

It's odd because the correct panorama shows if using Google directly. E.g.:

https://maps.google.com/maps?q=654+NW+4TH+ST,+REDMOND,+OR+97756

Here is a JSFiddle demonstrating the issue.


texts: Problem with KML in Google Earth iOS-App

texts:
I have a some what strange problem with KML/KMZ-files in the iOS google earth-app for ipad/iphone. I export my map from ArcMap to KMZ. Email myself and open the file with KMZ-loader and then open the file in Google Earth.

The problem is that google earth does not draw the points in the file. It works fine with lines and polygons and I can't understand what's the problem.

If I touch the "bookmarks" next to the search bar, I can see all the points and read the attributes for every point. I can even zoom to each one but nothing shows on the map.

Any one have a possible solution to this?

Regards/ Olof


texts: Polygon Extent to Xmin, Ymin, Xmax, Ymax .csv

texts:
I have a piece of code that reads in a .csv file and outputs a shapefile. After processing and altering the shapefile. The shapefile contains multiple polygons. I need to turn it back into the .csv, but am unsure how to "go backwards". I was thinking that I could do it somehow with the raster's resolution and polygons extents, but just am kind of confused with the process. Ultimately, I want to convert the .csv into an .xml.
Any advice?
training_csv = gpd.read_file('D:/data/farmer.csv')
training_csv.drop(['filename', 'width', 'height', 'class', 'geometry'], axis=1, inplace=True)

tif_file = 'D:/data/farmer/tif/0.tif'  # georeferenced TIFF
shp_path = 'D:/data/farmer/shape/'  # don't add extension

save_name = 0

# open georeferenced tif file
with rasterio.open(tif_file) as image:
    for i in training_csv.iterrows():
        xmin = float(i[1][3])
        ymin = float(i[1][2])
        xmax = float(i[1][1])
        ymax = float(i[1][0])

        # get vertices of the bounding box
        # geocoordinates from pixel coordinates
        p1 = image.xy(xmin, ymin)
        p2 = image.xy(xmax, ymin)
        p3 = image.xy(xmax, ymax)
        p4 = image.xy(xmin, ymax)
        print(p1)

        # save shapefile containing one bounding box shape
        w = shapefile.Writer(shp_path + str(save_name) + '.shp')
        w.field("name", "C")  # pyshp needs at least one field
        w.poly([[p1, p2, p3, p4]])  # generate bbox polygon
        w.record('bbox')
        w.close()

        # generate .PRJ file
        crs_wkt = image.crs.to_wkt()
        prj = open(shp_path + str(save_name) + '.prj', "w")
        prj.write(crs_wkt)
        prj.close()

        save_name = save_name + 1

Original CSV


texts: Optimize GEOSGeometry's spatial operations by avoiding thread-safety checks

texts:
A ranking algorithm I'm developing for a Django webapp requires performing spatial operations (within, intersections, sym_difference, etc.) using a query against large number of geometry objects.

Presently, the algorithm is very time-consuming, and the profiler reveals that the bottleneck is django.contrib.gis.geos.prototypes.threadsafe.__call__, which are called by virtually every methods in GEOSGeometry. From what I understand, threadsafe.py is a wrapper to the underlying GEOS's threadsafe API.

In this case the threadsafe-checking presents huge overhead to the computation time. Also my computation is strictly single threaded. From these reasons, I wonder is it possible to bypass or avoid the thread-safety checking performing by Django's GEOS API?


texts: Creating custom Editor Object for ArcMap and tie it to custom tool?

texts:
I have some spatial enabled tables tables in SQL Server 2008 and add them to ArcMap using custom layers(or a a spatial enabled OLE DB provider) successfully,but as you know,the editor toolbar of ArcMap doen't get enabled for these layers.
Is it possible (and if yes how) to create a custom "Editor Object" for ArcMap which implements all necessary interfaces(like IEditor,IEditor2,ISnapEnvironment,...) and use this custom editor for these custom layers?and is it possible to create a "Sketch Tool" which can work with this custom editor and have same properties like snapping?


texts: Creating rectangle over a buffer circle with specified angles/rotation?

texts:
I have a layer that is a circle buffer around a point. What I want to do is create a rectangle around the circle, but have every rectangle be perfectly lined up N/S (0degrees of angle/rotation). I have tried using Minimum Bounding Geometry, but the rectangles are all rotated a little different each time instead of all being the same. Is there any way to do this?

Thank you for the help.


texts: Only one country overlay at a time from table features in GEE

texts:
Can anyone help me to overlay the only country on the LST map on Google Earth Engine using the county features table?
// Import counties.
var counties = ee.FeatureCollection("TIGER/2018/Counties")
  .filter(ee.Filter.eq('STATEFP', '41'));
print(counties);

// Import LST.
var lst = ee.ImageCollection('MODIS/006/MOD11A1')
  .filterBounds(counties)
  .filterDate('2020-03-30', '2020-04-29')
  .select('LST_Day_1km');

// Add LST units as degrees c as new band.
var tempToDegrees = lst.map(function(image){
  var degreesC = image.multiply(0.02).subtract(273.15).rename('LST');
  return image.addBands(degreesC);
});

// Make a LST time series per county in FeatureCollection.
var lstTsByCounty = counties.map(function(county) {
  return tempToDegrees.map(function(img) {
    // Reduce precision of county boundaries to match the LST resolution;
    // avoids issues with exceeding memory limit.
    var geom = county.geometry().simplify(1000);
    var stat = img.select(['LST']).reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: geom,
      scale: 1000,
      bestEffort: true,
      maxPixels: 1e13,
      tileScale: 4
    });
    return ee.Feature(geom, stat)
      .copyProperties(county, ['NAME'])
      .set('system:time_start', img.get('system:time_start'));
  });
}).flatten();

// Export the LST time series per county FeatureCollection. The job is
// too big to handle interactively.
Export.table.toAsset({
  collection: lstTsByCounty,
  description: 'lstTsByCounty',
  assetId: 'gis_se_q_365561_68793'
});

// After the export finishes- load the FeatureCollection to view the chart.
// You'll need to change the path to your own account.
var lstTs = ee.FeatureCollection("users/nidhivermaiiita/gis_se_q_365561_68793");
print(lstTs.limit(5));

// Chart the LST time series of all the counties.
// If you don't like all of the counties in the same chart,
// filter the lstTs collection by the 'NAME' property.
var chart = ui.Chart.feature.groups({
  features: lstTs,
  xProperty: 'system:time_start',
  yProperty: 'LST',
  seriesProperty: 'NAME'
});
print(chart);

// Display LST to Map.
var landSurfaceTemperatureVis = {
  bands: ['LST'],
  min: -10,
  max: 35,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};
Map.centerObject(counties);
Map.addLayer(
  tempToDegrees.mean()
  .clipToCollection(counties), landSurfaceTemperatureVis, 'LST');

// Display the county boundaries.
var countyLines = ee.Image().byte().paint({
  featureCollection: counties,
  color: 1,
  width: 2
});
Map.addLayer(countyLines, {palette: '000000'}, 'countys');
Map.setOptions('SATELLITE');


texts:
I am trying to create a bubble map in CARTO in which a user can click on a bubble and one of the fields that appears is an mp3 file. Alternatively, a soundcloud page with the mp3 could appear with the bubbles' other fields.

My research on the subject has not been fruitful, although I have seen maps that do this successfully: http://www.antievictionmappingproject.net/narratives.html (here the blue bubbles have audio come up when they are clicked).


texts: cannot download 3rd party plugins in qgis

texts:
I want to download SDA4PP plug-in in QGIS. Unfortunately I cannot see the button to download 3rd party plugins, as indicated in the install instructions at that link. I am using this program on a company network. I have put in some proxy settings which mean that I can access the global/official repository. Before I did this I could not even access this. When I go to the plugins tab though there are only a few listed and when I put something in the search box nothing comes up, i.e. it only searches currently downloaded plugins (those that came with the install).

Anything I'm missing?


texts: Show geospatial TIFF as base layer in OpenLayers using GeoServer

texts:
I am trying to show a GeoSpatial TIFF file as a base layer in OpenLayers, on which further layers can be imposed.

Upon trying to add the TIFF file as GeoTIFF store GeoServer throws below error:

Could not list layers for this store, an error occurred retrieving them:

Failed to create reader from file:data/Sample.tif and hints Hints: REPOSITORY = org.geoserver.catalog.CatalogRepository@30c5d9f7 EXECUTOR_SERVICE = java.util.concurrent.ThreadPoolExecutor@45cb6ccc[Running, pool size = 0, active threads = 0, queued tasks = 0, completed tasks = 0] System defaults: FORCE_LONGITUDE_FIRST_AXIS_ORDER = true FORCE_AXIS_ORDER_HONORING = http FEATURE_FACTORY = org.geotools.feature.LenientFeatureFactoryImpl@3cc41abc LENIENT_DATUM_SHIFT = true GRID_COVERAGE_FACTORY = GridCoverageFactory TILE_ENCODING = null STYLE_FACTORY = StyleFactoryImpl COMPARISON_TOLERANCE = 1.0E-8 FILTER_FACTORY = FilterFactoryImpl

I tried to add the TIFF as ImageMossaic as well, but GeoServer gave below error:

Could not list layers for this store, an error occurred retrieving them:

Failed to create reader from file:data and hints Hints: REPOSITORY = org.geoserver.catalog.CatalogRepository@30c5d9f7 EXECUTOR_SERVICE = java.util.concurrent.ThreadPoolExecutor@45cb6ccc[Running, pool size = 0, active threads = 0, queued tasks = 0, completed tasks = 0] System defaults: FORCE_LONGITUDE_FIRST_AXIS_ORDER = true FORCE_AXIS_ORDER_HONORING = http FEATURE_FACTORY = org.geotools.feature.LenientFeatureFactoryImpl@3cc41abc LENIENT_DATUM_SHIFT = true GRID_COVERAGE_FACTORY = GridCoverageFactory TILE_ENCODING = null STYLE_FACTORY = StyleFactoryImpl COMPARISON_TOLERANCE = 1.0E-8 FILTER_FACTORY = FilterFactoryImpl

I have hundreds of continuous Geospatial TIFF images which shall be visible in OpenLayers according to the current position.

Any pointers/suggestions will help. Please suggest if I need to revise my approach or try any other way.


texts: Select all "NULL" values in the attribute table?

texts:
I have a quite simple syntax question.

I want to select ALL my null values in an attribute table of a feature class, not just the values in a single field.

So the question would be: what is the valid sql statement that could get me all the null values?

Thank you!


texts: Instalation GRASS add-ons and its usage in QGIS

texts:
How to properly install (any procedure!) and enable it to have the ability to run them via QGIS and use them for analyses in this environment?

There is lot useful algos for GRASS mentioned on the GRASS website but not available in QGIS. Let's say that tutorial is a bit outdated :)


texts: Calculating slope over distance from raster

texts:
I am sitting with an elevation raster, trying to produce a calculation for slope angle. I have a cell size of 2,-2 and I am trying to figure out how to calculate cell value difference between adjacent cells. Or is there another way to visualize the slopiest slopes in my raster?
I am currently using QGIS and SAGA.


texts: I need to use the Add Join tool but before I can do that I must get the csv file into ArcMap. Yet when I go to ArcCatalog to try to add the csv it is not there. Also I tried using simply the add data button but again the csv is absent. So how can i get the csv file to show up in ArcCatalog so I can begin the process with the Add Join tool?


texts: how to compute tile row from latitude in mercator projection

texts: I work on google tile images. I want to compute a tile row from its latitude. How can I do that. I know that tile system cover all of the earth from -85 latitude to 85 latitude. 4 tiles exists at the first level. Google uses mercator projection.


texts:
I have a file named 'VNP13A1...h5' and I'm trying to georeference it in Python 3. I know how to read and extract data from this file and what I'm trying to do is automating this process for other files as well.

This is what I did so far

from h5py import File
import rasterio as rio
import numpy as np

file    = File(i_file, 'r')
file_md = file['HDFEOS INFORMATION/StructMetadata.0'][()].split()
file_md = [md.decode('utf-8') for md in file_md]

for md in file_md:
    if 'UpperLeftPointMtrs=(' in md:
        mtrs   = md.split('=', 1)[1]
        ulc    = np.fromstring(mtrs[1:-1], dtype=float, sep=',')
    elif 'LowerRightMtrs=(' in md:
        mtrs   = md.split('=', 1)[1]
        lrc    = np.fromstring(mtrs[1:-1], dtype=float, sep=',')
    elif 'XDim=' in md:
        mtrs   = md.split('=', 1)[1]
        width  = float(mtrs)
    elif 'YDim=' in md:
        mtrs   = md.split('=', 1)[1]
        height = float(mtrs)

res   = np.divide(np.subtract(lrc, ulc), (width, height))
prj4  = '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs'
trans = rio.Affine.translation(*ulc) * rio.Affine.scale(*res)

For this specific file proj4 value is

+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs

I searched the metadata and found this

`file_md[:15]

['GROUP=SwathStructure',
 'END_GROUP=SwathStructure',
 'GROUP=GridStructure',
 'GROUP=GRID_1',
 'GridName="NPP_Grid_16Day_VI_500m"',
 'XDim=2400',
 'YDim=2400',
 'UpperLeftPointMtrs=(3335851.559000,4447802.078667)',
 'LowerRightMtrs=(4447802.078667,3335851.559000)',
 'Projection=HE5_GCTP_SNSOID',
 'ProjParams=(6371007.181000,0,0,0,0,0,0,0,0,0,0,0,0)',
 'SphereCode=-1',
 'GROUP=Dimension',
 'END_GROUP=Dimension',
 'GROUP=DataField']

What I can say about 'ProjParams=(6371007.181000,0,0,0,0,0,0,0,0,0,0,0,0)' is that I think 6371007.181000 is related to radius of the ellipsoid used in this projection but nothing else.

What are those other values used for?


texts:
Does anyone know where the Raster Clipper tool went in QGIS 2.18.20?

It used to be under Raster->Extraction->Clipper.


texts: Displaying features on JMapPane in GeoTools from postgresql database table geometry column?

texts:
I am making a GIS maven project in GeoTools.

How can we display data from the PostgreSQL database in GeoTools?

Presently I can display the shapefile from the file path, but I want to display this shapefile from the PostgreSQL database. So I could run spatial queries here in GeoTools. I gave all the dependency for connectivity of PostgreSQL database to GeoTools in pom.xml file.


texts: Is it possible to 'scrape' shapefile layer from a map saved as pdf/jpeg?

texts:
Data on the changes in voting district delimitations in SA is available as a series of maps (no shapefiles made official).
I'd like to create a shapefile by
scraping the boundaries from the individual pdf's and putting them together.
Is it viable?
What would be the easiest/fastest/best way?
Ideally the process would be automated.

Example of the document:
Schedule 1 MSA Detailed ward maps (last downloadable link)
http://www.demarcation.org.za/index.php/ward-delimitation/2014-2015/western-cape/city-of-cape-town


texts:
I have a big project (whole country) and it contains smaller features (administrative units).

My project is in .shp.

How to export features' coordinates one by one or all at once as long as all features' coordinates are separate / distinctive by eye in order to copy-paste them separately?

I don't care how coordinates are formatted as long as I can easily copy-paste them without downloading some huge program in order to open exported file.

But if there's any change to control the output, my desired final format is this:

//             point 1                           point 2
[["lat"=>xx.xxxxx,"lng"=>yy.yyyyy],["lat"=>xx.xxxxx,"lng"=>yy.yyyyy],/*etc*/]

Otherwise I could just program a JS function that formats the coordinates to my desired format.

All these features are on same layer:


texts: Does GeoServer support OpenStreetMap OSM format?

texts:
Apologies if the question has been asked but I am kind of lost between al the resources and tools available.

I set up a GeoServer and have been reading the documentation but I still don't understand how to serve tiles from OpenStreetMap.

I might have misunderstood the tools to achieve building a simple world map with multiple styles. Using Leaflet and OSM works fine as a demo app and as I understand, I need to have multiple maps generated on my server in order to have custom styled maps. The only layers I want are roads, boundaries, with different colors.

Am I using the wrong tools?
If not, how is GeoServer going to help me serve maps with the following example set of URLs and layers?

Example :

myserver.com/dark/z/x/y
myserver.com/blue/z/x/y
etc

At first I used Mapbox and created multiple maps but their pricing is too high for my small starting business, but it does what I want to achieve. I would like to do it myself and I have a developer background.

Again I'm sorry I'm mixing everything, I'm willing to get more documentation on the matter and I think this could be basic to achieve as I don't want to get deep into customization. Please feel free to redirect me to any tool/documentation in order to host my custom tile server.


texts: qgis-bin.exe - Entry Point Not Found (QGIS 2.4 / 2.5)

texts:
Since upgrading to QGIS 2.4 (Chugiak) I have been receiving the below pop-up error...

qgis-bin.exe - Entry Point Not Found
The procedure entry point xmlSchemaNewDocParserCtxt could not be
located in the dynamic link library C:\OSGEO4~1\bin\spatialite.dll

I have tried multiple re-installs, once with the bare minimum of software and libraries, but continue to get the same error. The required dll file IS in the bin folder and exists no where else on the system (I've read that should the dll exists ion the system32 folder it will act as the default and throw a wrench into the functioning of any dependant software).

I have also adapted the solution at https://gis.stackexchange.com/a/30450/12354, but to no avail.

Any advice?


texts: Creating Slope from DEM in ArcGIS Yields Unexpectedly High Values

texts: I have a DEM data set of Nepal and want to create a percent slope map of that data. When I use the slope function in ArcToolbox I am getting extremely high values (I know Nepal is very mountainous but i am sure that 500,000,000+ is incorrect).


texts: Selecting polygon features that have centers outside of another polygon feature class using ArcMap/ArcPy?

texts: Since I only have Incorporated area feature class, how can I select census blocks that are in unincorporated area? Or not within any incorporated areas using ArcMap or ArcPy?


texts: Blocking areas from consideration when running r.drain to calculate walking path from cost map?

texts:
r.drain was obviously originally intended to be used for drainage catchment, which is a problem when trying to find a least-cost walking path, as r.drain is happy to send you walking down a riverbed.

I've thought of blocking areas with either "nodata" or 9999 value cells, but i'm not really sure how exactly to do this.  Is there an easy way?


texts: Traces outside polygons after merging several polygons together

texts:
I'm using ArcMap 10.3.1

I have a polygon shapefile for areas upstream of 5 monitoring sites. I used dissolveand merge to merge all these polygons together. You can download the resulting shapefile from here.

I used this shapefile in R v. 3.2.5 for further analysis. In the result of this analysis, below, there was a line cutting through the polygons

I checked the shapefile again in ArcMap 10.3.1 to find that after merging, there were traces outside the polygons.

I tried to remove them but I couldn't.


texts:
I have bought a few (30) GeoTiff maps from my country's Geo-department - South Africa. They all have a date in the filename and the newest one was created in 2004, I see one is dated 1977 - that old! Anyway, they came with a .tfw file each. When I simply copy these into Google Earth, or use MapTiler to create a default-settings WSG84 tile set, the map is unaligned. Meaning that they are really close to their correct positions but a bit off center and the deviation from supposed-position differs everywhere on the image/map.

Either the maps are wrongly geo-tagged or I use the wrong formats. Please note: I have no clue or experience with maps and their types (WSG84/UTM etc) I need to convert these maps into tiles for my iPhone app.

I have also tried using maps of USA and they worked and were outputted perfectly.

All I can tell about these maps which may help pro's is:

It is a big white block with the map inside and the key/legend on the sides/bottom.
The map inside has an arced top and bottom and sides/longitude-lines are skew.

I would like to think it is the same as this image: (Only a small part thereof of coarse!)

I have tried to realign/resize the map as an image-overlay on Google Earth, but then had no idea of how to output it or get the needed (and corrected) data from it.

If anybody know how to use these, convert them to tiles or tips/tutorials that could help me fix this and get the maps straight, please tell!

If possible I would prefer to do it in MapTiler or at least on a Windows OS machine

Thank you!


texts: How to interpolate a raster surface from linear vector input?

texts:
Using ArcGIS 10.4.1 I have aircraft flight line data in vector format, with attributes for altitude, speed and rate of climb.  From those attributes, I have  estimated a relative loudness for each arc segment, stored in a new attribute field.

Now I need to estimate that loudness on the ground.  I do not need to calculate an exact decibel level; rather a simple, relative intensity will suffice.  I will assume that the sound decays over distance following the inverse square law, independent of frequency. I envision using raster analysis, where the output cell values derive from the arc segment loudness, decreasing in distance from the flight lines.

While there are various surface interpolation tools available (IDW, Kriging, etc), they all rely on a point input, whereas my input data is linear.  Is there an analysis approach that considers a linear source?


texts: ArcMap 10.3 Mosaic to New Raster produce a single color output

texts:
Here's what my final result looks like

The workflow is quite simple but I cannot get it go right

Create first the mosaic dataset

arcpy.CreateMosaicDataset_management(in_workspace="<path>/xxx.gdb", in_mosaicdataset_name="B01", coordinate_system="PROJCS['WGS_1984_UTM_Zone_31N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['false_easting',500000.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',3.0],PARAMETER['scale_factor',0.9996],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]];-5120900 -9998100 10000;-100000 10000;-100000 10000;0,001;0,001;0,001;IsHighPrecision", num_bands="", pixel_type="", product_definition="NONE", product_band_definitions="")

I then add my two adjacent (sslightly overlapping) rasters

arcpy.AddRastersToMosaicDataset_management(in_mosaic_dataset="<path>/xxx.gdb/B01", raster_type="Raster Dataset", input_path="B01_01;B01_04;",update_cellsize_ranges="UPDATE_CELL_SIZES", update_boundary="UPDATE_BOUNDARY", update_overviews="NO_OVERVIEWS", maximum_pyramid_levels="", maximum_cell_size="0", minimum_dimension="1500", spatial_reference="", filter="#", sub_folder="SUBFOLDERS", duplicate_items_action="ALLOW_DUPLICATES", build_pyramids="NO_PYRAMIDS", calculate_statistics="NO_STATISTICS", build_thumbnails="NO_THUMBNAILS", operation_description="#", force_spatial_reference="NO_FORCE_SPATIAL_REFERENCE")

Statistics:

arcpy.CalculateStatistics_management(in_raster_dataset="<path>/xxx.gdb/B01")

Color Balancing (if I do not do this, I got an end result where there's no unique value, but the color difference in both image is too bad)

arcpy.ColorBalanceMosaicDataset_management(in_mosaic_dataset="<path>/xxx.gdb/B01", balancing_method="DODGING", color_surface_type="FIRST_ORDER", target_raster="", exclude_raster="", stretch_type="NONE", gamma="1", block_field="")

arcpy.BuildOverviews_management(in_mosaic_dataset="<path>/xxx.gdb/B01", where_clause="", define_missing_tiles="DEFINE_MISSING_TILES", generate_overviews="GENERATE_OVERVIEWS", generate_missing_images="GENERATE_MISSING_IMAGES", regenerate_stale_images="REGENERATE_STALE_IMAGES")

The part where I create the unique output end document

arcpy.MosaicToNewRaster_management(input_rasters="<path>/xxx.gdb/B01", output_location="<path>/xxx.gdb", raster_dataset_name_with_extension="xxx_B01", coordinate_system_for_the_raster="", pixel_type="16_BIT_UNSIGNED", cellsize="", number_of_bands="1", mosaic_method="LAST", mosaic_colormap_mode="FIRST")

Creating Statistics is of no help

arcpy.CalculateStatistics_management(in_raster_dataset="<path>/xxx.gdb/xxx_B01")


texts: Delete tags from planet.osm with osm2pgsql and minimize size of database

texts:
I want to import into PostGIS the whole planet but I am interested for certain tags. Specifically I am not interested in nodes, ways and areas which are related with the land. Below is a small part from the openstreemap-carto.style

# OsmType  Tag          DataType     Flags
node,way   access       text         linear
node,way   addr:housename      text  linear
node,way   addr:housenumber    text  linear
node,way   addr:interpolation  text  linear
node,way   admin_level  text         linear
node,way   aerialway    text         linear
node,way   aeroway      text         polygon
node,way   amenity      text         polygon
node,way   area         text         # hard coded support for area=1/yes =>      polygon is in osm2pgsql
node,way   barrier      text         linear
node,way   bicycle      text
node,way   brand        text         linear
node,way   bridge       text         linear
node,way   boundary     text         linear
node,way   building     text         polygon

After that there is some information regarding the flag delete

 # Deleted tags
 # These are tags that are generally regarded as useless for most rendering.
 # Most of them are from imports or intended as internal information for mappers
 # Some of them are automatically deleted by editors.
 # If you want some of them, perhaps for a debugging layer, just delete the lines.

My question: If I use the flag delete for some nodes, does it mean that these nodes (or ways) will be completely discarded resulting in a smaller database?
Do I also need less space during importing planet in PostGIS, if I use the delete flag?

If not (I mean if the space needed during the process of importing data into db is the same and only the resulting database is smaller) suggest me a tool for dropping some nodes or ways with certain tags from the planet.osm in order to produce a smaller planet file for use with osm2pgsql.


texts:
Say I have the coordinates:
N: 5413457.7
E: 1528677.3

If I were to add 10 to the northing value would that make the new coordinates point 10 meters north of the original point?


texts: Extracting numbers and text from string field using Python Parser of ArcMap Field Calculator?

texts:
This is a similar question to Extracting numbers from string field using Python Parser of ArcMap Field Calculator?
I have a similar issue in ArcMap 10.6. I have a field named 'SIZE' that holds number values and size unit (e.g. '200 mm' or '12"').
I am trying first to extract the numbers and deal with the units after.
I created a new string field named 'SIZE2', while using the same code, but when I run it the values in 'SIZE2' stay NULL.
def makestr(SIZE):
 numlist = []
 for s in SIZE:
     if s.isdigit():
         numlist.append(s)
 return ''.join(numlist)

Any ideas?


texts: Converting from DMS to DD using Python in Field Calculator?

texts:
I need to convert lat/long that is expressed as degrees, minutes, and seconds in the data into decimal degrees. For example, in the data, they are listed as N335042.06 in the Latitude column, and W86031.04 in the Longitude column. I have done this problem before where I created a script which converted DMS to DD, and vice-versa, so I guess I could use bits from that. But the problem I'm having is how to ignore (for lack of a better word) the 'N' and 'W' in the data? Can I skip them? And DMS are listed all together without any symbols or spaces.

Can I use len(), range(), split() to specify what part to read from the value? For example, can do the following?

N335042.06 where,
33 = degrees
50 = minutes
42.06 = seconds ...?

I came across this ESRI article, but it's in VB. Will probably use it as a reference but some of the terminology/syntax is different from Python.

Final code that works!

# Pre-logic
def latDD(x):
  D = int(x[1:3])
  M = int(x[3:5])
  S = float(x[5:])
  DD = D + float(M)/60 + float(S)/3600
  return DD

# Expression
latDD(!Latitude!)


texts: Using Function with replace with the US-350 to 350 How to?

texts:
I am trying to replace the US- 350 with 350 in the Label Expression but I am stuck how to do that...

I am looking at the help here for the function:
Function FindLabel ( [LABELFIELD] )
FindLabel = replace([LABELFIELD], ", ", vbnewline)
End Function

This is what I did this :

Function FindLabel ( [ID] )
FindLabel = replace ( [ID], US- 350, 350)
End Function

no luck ??


texts: Postgis, get the points that are x meters near another point, in meters

texts:
I'm trying to find the points within a certain distance, in meters, from another point.

I'm using the function ST_Buffer, but I don't get how to express the distance in meters.

I'm trying with something like this:

select
1 as cartodb_id,
ST_Buffer(
  ST_Transform(
    ST_GeomFromText(
      POINT(-58.38145 -34.60368)'
      , 4326
    )
    , 3857
  )
  , 500
) as the_geom_webmercator

I have the coordinates from google maps (in wgs84, 4326)
then I transform them to webmercator (3857)
and then I pass 500 as parameter, and I get something that looks like the desired area, but after comparing it with what google maps says it falls short for a non neglectable distance.

Is this the right way to achieve it?

note: you can try the sql query in cartodb

--

From this question: https://gis.stackexchange.com/a/44481/19461 I came with the following code (I use SRID 31997 to calculate the buffer, and then go back to webmercator)

select
1 as cartodb_id,
ST_Transform(
  ST_Buffer(
    ST_Transform(
      ST_GeomFromText(
        'POINT(-58.427185 -34.624592)'
        , 4326
      )
      , 31997
    )
    , 2000
  ), 3857
) as the_geom_webmercator

Now if falls short for no more than 20 meters. I thought mercator (EPSG 3857) was in metric units.


texts: Projection of delimited text file does not match shapefile

texts:
I am using Windows 10 with QGIS 2.10.1 Pisa, and have tried the same on other QGIS versions with the same result.

I am trying to project a delimited text file onto another shapefile, however the two projections will not match whatever I do.

So far I have tried:

Using a background map (ne_10m_admin_0_countries) with WGS84:EPSG4326
Background map with WGS84:EPSG3857
Background map with conic map EPSG:102022
I have tried saving the text file as a shapefile in each case
I have also tried all of the above with the openlayers plugin
I have tried setting the same CRS of the shapefile to the text/shapefile

In all cases, the coordinates move around, but never fit to where they are supposed to. The other shapefile fits beautifully in each situation, and the GPS points I am trying to project have been checked and are correct.

Plot_code   Longitude   Latitude    Alt Al-S    As-S    Ba-S    Ca-S    Co-S    Cr-S    Cu-S    Fe-S    Hg-S    K-S Mg-S    Mn-S    Mo-S    Na-S    Nb-S    Ni-S    P-S Pb-S    Rb-S    S-S Sb-S    Se-S    Si-S    Sr-S    Ti-S    U-S V-S Y-S Zn-S    Zr-S
2-1-01-03-CT    0.53237 34.21674    1303.6  67370.08    3.78    375.26  8808.76 7.25    32.07   16.36   18019.95    307.41  48147.05    1832.97 681.62  2.19    7453.65 47.22   14.66   962.82  25.91   391.9   482.89  4.5 2.32    305318.2    143.92  2527.65 17.54   52.98   23.37   122.84  396.83
2-1-01-03-FT    0.53192 34.21778    1295.7  74421.06    4.28    299.04  6073.86 7.85    25.85   16.77   21118.91    2.54    47335.6 491.87  726.28  1.98    7225.7  52.78   18.44   564.39  29.49   417.32  178.11  4.6 2.15    282504.8    116.83  2624.23 19.16   57.39   26.66   78.01   464.67
2-1-01-03-MT    0.53234 34.21722    1300.1  72070.97    5.15    437.69  6356.46 7.63    45.29   12.44   19592.15    NA  51437.09    1187.14 697.89  2.18    7497.35 51.43   16.85   655.4   29.7    425.11  224.67  4.66    2.35    313366.77   113.77  2771.96 21.62   58.22   26.31   72.76   449.88
2-1-01-06-CT    0.53234 34.21628    1302.6  70188.74    3.92    309.57  6229.26 7.71    50.31   15.02   18467.58    NA  54242.03    481.89  477.39  2.12    7445.71 48.51   14.28   518.77  25.58   429.22  236.01  4.77    2.02    308015.37   110.12  2507.86 18.47   56.46   22.39   76.5    397.82
2-1-01-06-MT    0.53244 34.21593    1303.6  75760.14    5.17    420.7   5412.34 8.15    27.1    14.28   21409.93    NA  50840.38    89.47   548.52  1.99    7293.96 47.83   14.89   433.44  34.1    458.1   210.27  4.64    2.7 285652.42   120.56  2666.82 22.17   61.04   26.16   66.54   476.81
2-1-01-08-CT    0.53721 34.21689    1338.8  67282.55    4.39    340.27  5053.28 7.25    31.92   10.71   17658.11    NA  56488.27    1139.49 571.27  1.94    7347.79 52.95   13.45   339.38  28.37   436.39  130.16  4.68    2.36    325334.21   90.64   2617.23 19.46   56.02   24.8    57.74   439.38
2-1-01-08-MT    0.53656 34.21752    1328    63928.37    4.05    435.14  5353.69 6.9 41.04   10.33   15110.59    NA  56868.25    1137.7  544.17  1.99    7379.4  43.2    13.54   334.23  27.41   431.99  137.67  4.6 1.95    346056.7    91.98   2301.29 20.65   50.06   21.07   52.77   397.91
2-1-04-07-CT    0.59156 34.24355    1179.5  6901.13 2.2 436.73  2497.35 5.02    23.46   11.96   5719.69 NA  27490.48    3331.05 141.24  2.43    6654.78 30.8    5.38    NA  13.38   118.87  NA  4.63    1.78    548831.16   77.48   1550.86 4.87    39.25   10.19   28.88   724.74
2-1-04-07-FT    0.59153 34.24282    1179    487.21  2.01    377.85  2351.34 4.62    18.8    4.96    4649.69 NA  28050.41    3265.53 174.52  2.65    6497.28 24.99   6.38    NA  13.71   114.74  NA  4.82    1.31    565137.02   77.02   1477.41 5.33    38.45   13.59   26.12   653.12
2-1-04-07-MT    0.59153 34.24294    1179.1  8996.03 2.18    468.3   2580.06 5.04    24.49   8.77    5977.46 NA  28001.39    3198.4  164.27  2.48    6475.02 30.14   6.96    NA  13.4    117.33  NA  4.87    1.29    540065.51   75.52   1595.78 3.65    37.14   14.87   29.26   723.53


texts: How does integrate work with point data sets?

texts:
I'm trying to combine point data into clusters and I think I've come up with a method that will work for what I want to accomplish.  However, I don't understand why it works the way it does and I'm hoping someone can explain it to me.

Basically I have a set of wells that are operated by different companies and are scattered throughout PA.  I want to determine how clustered the wells are by operator in terms of how many clusters there are for each operator and how many wells are in each cluster (clustering at a distance of 5 miles).  It appears that if I run integrate on my point data for the wells it will do this.  Then I can collect event values to determine how many wells are in each cluster.  What I don't understand is why some of the data is broken into more than one cluster.

[Image Description] The three highlighted points are the integrated points for the company represented by the light yellow well sites.  When you open the attribute table for the integrated data it is clear that the group on the right is one cluster and that integrate has determined that the group on the left is two clusters split on the blue line.

What I don't understand is why the points immediately above and below the line have been assigned the clusters that they have been.  They are all within 5 miles (the integrate thresh hold) of the lower integrate point.  And the upper ones are more than 5 miles from the upper integrate point.  So why split it there?

Any thoughts would be helpful.

To help clarify I've added a screenshot of the process that I used to integrate all the points into the clusters that I would like to use.

http://imgur.com/swzxW


texts: Python - Getting coordinates without restrictions

texts:
I'm using Python for getting coordinates (latitude & longitude) for a comprehensive list of cities and postal codes (thousands of rows, nearing 10 000 rows).
At the moment I'm using geopy and Nominatim for downloading the coordinates, but Nominatim has restrictions for downloading large datasets.
My script is now constantly failing and I expect it to be because my own dataset of the addresses and postal codes is too large.
Is there a (free?) service, to be used with Python, that does not have a restriction on how many requests I am allowed to do with it / how many requests per second it allows? Google Maps or similar maybe?


texts: The WMS Access URLs and Data Format for viz.globe.gov

texts:
As I'm reading there are two popular formats for GIS data: KML and WMS. This tutorial says the access URL for the GLOBE WMS GetCapabilities is at http://viz.globe.gov/viz-bin/wmt.cgi?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetCapabilities. However the /viz-bin/wmt.cgi does not exist anymore.

The questions:

Has the WMS URL moved?
Has the data format changed? If it does,what is the new format of this site?

Edit: More questions:

Is the site still good for testing WMS service against?
What is the other best server for testing?


texts: generating some random points along network

texts: I'm having  trouble with generating some random points along a network. I am trying to do it with Arcobjects in the VBA environment. I would be grateful if you could help me.


texts:
Please help me on how can get my map from Geoserver. I want create my country web but I can nott see it in my browser but I can see in Geoserver. Please help me I am very new openlayers and Geoserver.

  <title>OpenLayers Example</title>


texts: How to randomly subset X% of selected points?

texts: What methods are available in ArcGIS 10.2 to randomly subset a selection of points.  For example, in the attached screenshot I am interested in keeping 20% of the selected points and deleting the rest.


texts: Cell size in Mosaic to New Raster versus Resample cell size

texts:
I have several TIF 2 meter LIDAR DEM files (100+) that I am merging into a single raster and then creating 4m, 6m, 10m, and 30m DEM raster variants.

I would normally just use the Mosaic to New Raster tool and not specify a cell size, maintaining 2 meters and then use the resample tool to downgrade the 2 meter DEM to the variants of 2,4,6,10,and 30 meters.

There is a Cell Size (optional) within the Mosaic to New Raster tool. If I input my 100+ 2 meter DEM files and then select a cell size of 4,6,10... etc., does this tool resample the DEM?

I am not clear on what this field does if I am creating a lower resolution raster from a higher resolution. If it does resample, what method is used since I cannot specify?

It seems I am safer now just Resampling. I would love to learn/understand more and save myself a step.

I have researched this, but I come across questions and answers around merging DEMS but not specifically addressing my question.

I am using ArcMap 10.1. There are 100+ TIF DEM files that are 2 meter cell size. 32 bit float and 1 band. Not sure if any of that is needed.


texts: polishing qgis map for publication using inkscape

texts:
I have a problem modifying the map I created in QGIS. I need to add some text on the map and axis labels. Apparently, labels cannot be rotated in qgis, so this I need to do in inkscape or in a related programme.

If I export the map as pdf, I cannot adjust the font size for the "tick marks" (lat and long values) according to the guidelines of the journal I want to submit my article to.
I discovered the simplesvg tool, which gives a very nice resolution map, but then I don't have the coordinates along the borders of the map (or is there a way with this tool?)

Is there someone that can help me? I'm getting pretty desperate.. :s

Thanks!
Ellen


texts: QGis Standalone App access to maptools defined in app directory

texts:
I'm developing an QGIS application estandalone in c ++, and would like to know why the "app" directory, the QGIS source code is not exposed. As I would like to use some of the Tools that there are defined.
It can be done in any way? Which would be the most elegant way to do this?


texts: SLD specification for a vertical line on every line turning point

texts:
You can define a vertical line in SLD like this:

   <LineSymbolizer>
     <Stroke>
       <GraphicStroke>
         <Graphic>
           <Mark>
             <WellKnownName>shape://vertline</WellKnownName>
             <Stroke>
               <CssParameter name="stroke">#333333</CssParameter>
               <CssParameter name="stroke-width">1</CssParameter>
             </Stroke>
           </Mark>
           <Size>12</Size>
         </Graphic>
       </GraphicStroke>
     </Stroke>
   </LineSymbolizer>

There is no definition in xml for the place of that vertical line. I am making a solution for OpenLayers like in [1] but cannot figure out, how to place the vertical lines across the long line that goes horizontally.

Adding to OpenLayers goes similar to:

var format = new OpenLayers.Format.SLD();
OpenLayers.Request.GET({
    url: "yoursldfile.sld",
    success: sldparser
});
function sldparser(req) {
    sld = format.read(req.responseXML || req.responseText);
    styles = sld.namedLayers.interpreted.userStyles[0];
    building_vec.styleMap.styles.default = styles;
}

I would like to add the vertical lines to places, where the continuous line makes a tilt to new direction.

Question:

How to tell the place of vertical lines in the combo of Openlayers and sld?

More specifically:

How to tell it is in the turning point of straight line?

With picture:

I want the vertical "railroad hatches" to appear only in places, where the "railroad" changes direction. [2]

My source:

[1] How to use a sld file with Openlayers using OSM data?

[2] http://docs.geoserver.org/stable/en/user/styling/sld-cookbook/lines.html#railroad-hatching


texts: Assigning coordinates to string in python?

texts:
Is there a specific Python package that allows me to assign coordinate location to a string?

For instance, say I have location data in the form of names, such as Burbank, CA.

Is there a module that allows me to assign coordinates to that string?

I'm new to programming so I'm not aware of all its functionality.


texts: QGIS OGR error creating failure - failed to prepare SQL

texts:
I'm working in a QGIS v3.4.8 Project with a number of layers including vectors in .tab and .gpkg formats and a raster aerial image. I'm editing .gpkg layers by creating either points or polylines. The creation and saving works fine for quite a while but then I am unable to save the layer/s.

I get a

  WARNING Commit errors

followed by a

  CRITICAL OGR error creating failure : failed to prepare SQL

QGIS is sometimes then not responding.

I'm on a Windows 10 machine. Sometimes I can toggle off the editing and then toggle on and the save will work. On other occasions I have to close down QGIS as it remains unresponsive. Does anyone have an idea why this is happening?


texts:
I just installed Windows 10, and I have ArcGIS 10.3. Now, when I use any script to create a new .mxd map document (and occasionally if I download an .mxd made on another computer), it shows up as a Microsoft word file, i.e. MAP.mxd.doc. I can't open it with ArcMap (even by searching for another program to use it, and selecting ArcMap as the program).

Does anyone know how to deal with this problem?


texts: QGIS - mbtiles export in core functionality?

texts: Seems like exporting mbtiles is a fairly common need, and it seems like there is a patchwork of various third-party pieces of code and old or partially-integrated plugins to address the need with various levels of success and quality - but they are all obscure and difficult and not-quite-there despite the developer's best efforts...  has the idea of building mbtiles export functionality into the qgis core ever come up?


texts:
Of recent I have this problem with Qgis 2.14 that my map composer won't produce a map of my project but instead it goes all blank. I even used 2.12 and even after uninstall/install the problem still persists.
What could possibly the problem?
I run a windows 10 OS on an i3 hp


texts: How to re-use an expression, stored as a string, in QGIS Field Calculator

texts:
As discussed in this question, QGIS automatically saves recently used expressions in the field calculator. But recently used expressions are deleted after a while. Sometimes I need to re-use an expression that I used a while ago, and it's no longer available.

So I needed a way to store expressions for re-use in the Field Calculator at a much later time.

The method I figured out is to store the expression as in a text field. When I need to re-use the expression, I copy it from the attribute table into the Field Calculator. Here's what that looks like with a simple example:

Of course it's trivial to remember and re-enter the expression in this example. I would use this method with much more complicated expressions, such as "Population"/$area*3.58701e-8/"Population", which calculates population density (people per square mile) after first converting square feet to square miles.

Now I'm looking for a better way than manually copying the expression. It seems like there should be a way in the Field Calculator to take the field name as input, interpret it as an expression, and output the result of that expression. Here's what I've tried so far:

Just using the field name outputs the formula as a string.

Is there a function that will force the Field Calculator to interpret a string as an expression? By comparison to the expressions to_real() and as_string(), I'm thinking it would be something like as_expression() or as_formula(). So the output of as_expression('$area') would be the calculated area of the polygon feature.


texts: How can I stop 'drag' moving my whole map in Leaflet.

texts:
I have a leaflet map, containing a tileLayer, and a layerGroup

In the layerGroup, I plot circleMarkers, and then connect them with a polyline.

Though the markers are defined as draggable, I get no drag events for them, and any attempt to drag them results in moving the whole map. I thought it must be because the polyline is on top of the markers, but no such luck. polyline.bringToBack() didn't help.

var multiLine = [];
var line = [];
towTrack.clearLayers();
for (var ix = 0; ix < display_rows.length; ix++) {
    var row = display_rows[ix];
    var lon = row.EventLongitude;
    var lat = row.EventLatitude;
    if (lat != null && lon != null) {
        var pos = [lat, lon];
        line.push(pos);
        var marker = L.circleMarker(pos, {
            radius: 5,
            color: 'yellow',
            weight: 1,
            draggable: true
        }).addTo(towTrack);
        marker
            .on('dragend', function (e) {
                var position = e.target.getLatLng();
                e.target.setLatLng(new L.LatLng(position.lat, position.lng), { draggable: 'true' });
            });
        if (row.EventCode == 4 || row.EventCode == 9) {
            multiLine.push(line);
        }
    }
}
var polyline = L.polyline(multiLine, { color: 'red', weight: 1 }).addTo(towTrack);
polyline.bringToBack();


texts:
It looks like the only way to get flat files/data from a WPS service is to use a client.  I'm looking to create a super-simple WPS client in Python 2.6 (or other language) without any dependencies on non-default modules.

Right now I'm gathering resources - I haven't done a lick of coding yet.

Can I use a POST HTTP request with a proper header to request a WPS file?

Additional Info:

There's an OWS LIB example client, but that requires an OWSlibrary.
There's a Javascript/OpenLayers client from PyWPS  - this looks like its the most independent of external dependencies sans the openlayers.js, but then ... well, OpenLayers is its own client.
I have the Firefox HTTP Resource Test for POST and GET requests
I have FireBug


texts: MapProxy installation on Windows: "SyntaxError: Non-ASCII character '\x90' in file"

texts:
I having a problem installing MapProxy 1.9.0 on Windows 7. I have Python 2.7.12 installed. I followed the instructions from this documentation.

But when I try to launch MapProxy, I get the following error:

(mapproxy_venv) C:\mapproxy_venv\Scripts\mymapproxy>mapproxy-util serve-develop mapproxy.yaml
[2016-11-10 09:47:57,851] mapproxy.config - INFO - reading: C:\mapproxy_venv\Scripts\mymapproxy\mapproxy.yaml
[info]  * Running on http://127.0.0.1:8080/ (Press CTRL+C to quit)
[info]  * Restarting with reloader
  File "c:\mapproxy_venv\Scripts\mapproxy-util.exe", line 1
SyntaxError: Non-ASCII character '\x90' in file c:\mapproxy_venv\Scripts\mapproxy-util.exe on line 1, but no encoding declared; see http://python.org/dev/peps/pep-0263/ for details

mapproxy-util --version: MapProxy 1.9.0

python -V: Python 2.7.12


texts: QGIS point sampling tool error in extracting values

texts:
I am trying to use the point sampling tool plugin from QGIS to get some sea ice data (ranging from 0-1) from a GTiff raster layer (ESPG: 3411) for my sites (ESPG:4326) in the Arctic. But the sea ice values obtained are all 1.0000.

I have set the project CRS as ESPG: 3411 (with OTF), same as the raster layer.

The overlay looks fine as I can actually click individual points on the QGIS to obtain sea ice data value point by point but I have hundreds of points.


texts: How to Select features greater a given size

texts:
My question is simple, but I didn´t find any tool in ArcMap that is able to achieve it. I have a quite large dataset whose features I want to work with. In order to handle the limited capacities of my computer I split the entire working-area (that is, the envelope of all my features) into a grid of size 5km. The idea behind this is that whenever I finish processing one of those grid-cells I can completely release their memory. However my grid-cells still seem to be to big as I still have memory-issues.

So I further thought that I could reduce the grid-size to the size of the biggest feature, that is either the one with the biggest width or height. The area however isn't useful to me, as there are also polygons that are very long but have a small area (street-features for instance).

The problem now is: How do I determine the largest feature to have an indicator for my grid-size? I found some arcpy that extracts the features envelopes:

arcpy.MinimumBoundingGeometry_management("inFeatures", "Envelopes.shp")

But now I´m stuck.


texts: I have a JPEG RGB map image composed of 12 colours produced by a programme using a colour palette. It converts correctly to  a modified TIFF, which prints out correctly. But when this is loaded into QGIS it appears in the canvas and print-out with the two green shades rendered as white. How can I correct this?


texts: C# Library for converting Shapefile Projection

texts: Is there any C# library which can take a shapefile as an input with projection information defined (e.g, EPSG 3011) and return a new shapefile with user specified projection (e.g, EPSG 4326) ?


texts:
I currently have a layer in GeoServer which consists of a number of PostGIS geometries. They are given color by a color column which is read in as part of the GeoServer style.

I now need to add tooltips to some of the shapes, preferably by adding a new column to the PostGIS table with the geometries and colors. Is there a way to do this using only PostGIS and GeoServer?

As a simplistic example assume I had a series of dots representing major cities, and I wanted a user to be able to move a mouse over the dot in openlayers to see the "City, Country" in a tooltip. In some cases that information may not be available (is null), in which case no tooltip should be shown.

I was thinking the answer was to create a TextSymbolizer label, but don't know how to get that to hide until mouseover.


texts: GeoPandas does not identify points inside polygon

texts:
I am trying to intersect some points with a polygon in GeoPandas.

First I create a Point GeoDataFrame from a text file.

import geopandas as gp
import numpy as np
import pandas as pd
from shapely.geometry import Point, Polygon

tile_index = pd.read_csv('../tile_index.dat', names=['tile', 'x', 'y'], sep=' ')
geometry = [Point(r.x, r.y) for ix, r in tile_index.iterrows()]
tile_index = gp.GeoDataFrame(tile_index, geometry=geometry)
tile_index = tile_index[(tile_index.x.between(-10, 30)) & (tile_index.y.between(-10, 30))]
print tile_index.head()

Which produces:

     tile      x       y                                        geometry
202   226 -9.226  -0.724  POINT (-9.226000000000001 -0.7240000000000001)
203   227 -9.226   9.276                POINT (-9.226000000000001 9.276)
204   228 -9.226  19.276               POINT (-9.226000000000001 19.276)
205   229 -9.226  29.276               POINT (-9.226000000000001 29.276)
219   243  0.774  -0.724               POINT (0.774 -0.7240000000000001)

Then I create a Polygon to intersect:

poly = Polygon([(0, 0), (0, 20), (20, 20), (20, 0)])
P = gp.GeoDataFrame(data=[1, poly]).T
P.columns = ['ID', 'geometry']
print np.any(tile_index.within(P))

Which returns False but at least 4 points should be within.

If I create a point using the coordinates of point 237 it does intersect.

pt = Point(10.774, 9.276)
p = gp.GeoDataFrame([[1, pt]])
p.columns = ['ID', 'geometry']
print p.within(P)

Which returns True.

If I do P.contains(tile_index.geometry) then all are False but if I do np.any([P.contains(row.geometry)[0] for ix, row in tile_index.iterrows()]) then this returns True.


texts: Cellular Automata Markov for forecast LULC in R software

texts:
I would like to forecast land use changes by implementing the Cellular Automata (CA) Markov Chain models in R software.
I have the classified LULC maps from 3 times periods (1992,2003,2014), and also the Markov transition probabilities matrices. However, I  was not able to find any implementation in R that project land use using the CA Markov chain. I would like to predict for 2025.

Does anyone know if is any implementation in R for projecting LULC, similar to the one implemented in IDRISI with the module (CA-Markov)?

Or anyone knows a way of predict LULC in R using the CA-Markov?


texts: Unable to open shapefile in ArcMap. "Error: Failed to open the dataset." Question mark icon and "ml" prefix on file name

texts:
ArcMap and ArcCatalog both show the a "?" next to several shapefiles in the catalog window. Additionally, a two-letter prefix "ml" has been added to the file name, but only when viewed in Catalog, not in Windows File Explorer. I am unable to open, load or even view properties of these files. Coworkers do not have this problem with the same files (in a networked shared folder). There are other files in the same folder that do not have this issue. This question is related, but the problem there was a missing .dbf file. All files that I have showing this error contain files with extensions: .cpg, .dbf, .prj, .sbn, .sbx, .shp, .shp.xml, .shx, so this didn't solve my issue.

What I have tried:

Adding file to ArcMap .mxd. Got error "Error: failed to open the dataset"
Opening file in ArcCatalog. Got error "Error: failed to open the dataset"
I was able to successfully open and view the same file from within Python/Geopandas with no problems, so I assume this is an ArcGIS-specific issue.
I can open and use any of these shapefiles in QGIS without any issues.

I'm running ArcGIS version 10.5.1 on Windows 10.
What am I missing?
EDIT based on comments:
The folder where these files are located, "County Layers (check against Shared Data)", is not a .gdb folder, nor is it in one. There are other .shp files in this same folder without this issue.
Here's all the component files of ParcelAddressOwner.shp as shown in File Explorer:

Based on Michael Stimson's comment, I used Geopandas to write out a new .shp file. I wrote the file, titled "Zone A Project Creeks - wtrcrs_line_Test.shp", to 3 different folders: 1) the original folder "County Layers" 2) one level up in "Data" and 3) two levels up in "GIS". Note that the file saved to folders "Data" and "GIS" in the image below do not have the ? or "ml" prefix.
Here's the full catalog window of these folders and test files.

This suggests a problem with the folder itself, but I'm not sure what. It has the same permissions as other folders at this level and above it.


texts: Why is overview with only one level larger than original GeoTIF

texts:
I have a large GeoTiff that I'm preparing for display through GeoServer. The original image has 4.9Gb and when I created an external overview with only one level (1/4 the size), using the same compression and color_space, the file size was 8.8Gb.

Shouldn't the overview be smaller?

Update: I ran the output from gdaladdo through gdal_translate, setting JPEG compression and YCbCr colorspace. File size went from 8.8Gb to 1.4Gb. Could it be that gdaladdo is not applying the compression? I have gdal v 1.11.3 - quite old. Will ask for our admin to update the package.

Here is the gdalinfo output of the original file:

Driver: GTiff/GeoTIFF
Files: sul2014_rapideye_b321_2sd.tif
       sul2014_rapideye_b321_2sd.tif.aux.xml
Size is 206279, 239176
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-57.754539576021727,-22.341123027138849)
Pixel Size = (0.000048007576766,-0.000048007576766)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  ( -57.7545396, -22.3411230) ( 57d45'16.34"W, 22d20'28.04"S)
Lower Left  ( -57.7545396, -33.8233832) ( 57d45'16.34"W, 33d49'24.18"S)
Upper Right ( -47.8515846, -22.3411230) ( 47d51' 5.70"W, 22d20'28.04"S)
Lower Right ( -47.8515846, -33.8233832) ( 47d51' 5.70"W, 33d49'24.18"S)
Center      ( -52.8030621, -28.0822531) ( 52d48'11.02"W, 28d 4'56.11"S)
Band 1 Block=1024x1024 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=1024x1024 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=1024x1024 Type=Byte, ColorInterp=Gray
  NoData Value=0

This is the gdaladdo command I used:

gdaladdo -r cubic -ro --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config BIGTIFF_OVERVIEW YES sul2014_rapideye_b321_2sd.tif 2

And here is the gdalinfo for the overview image

Driver: GTiff/GeoTIFF
Files: sul2014_rapideye_b321_2sd.tif.ovr
Size is 103140, 119588
Coordinate System is `'
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,119588.0)
Upper Right (103140.0,    0.0)
Lower Right (103140.0,119588.0)
Center      (51570.0,59794.0)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  NoData Value=0


texts: Zonal Statistics - polygons and rasters with different projections

texts:
I know that there are many many posts about this type of topic already, but all solutions did not work out for me so far.
I have a shapefile with polygons, which is described here and can be downloaded here (there is also a paper published in the Journal of Urban Economics, see this link). It comes as a geopackage. As the first link shows, the projection is World Mollweide (EPSG: 54009). I want to calculate zonal statistics for these polygons with various raster files, e.g. air quality (downloadable here (as .nc) or here (as .asc), which use the WGS84 projection. However, I just cannot manage to align the raster files with the shapefile.
I tried to save the polygons as a shapefile, assigning to it the WGS84 (EPSG: 4326) CRS. However, using QGIS 3.8 this results in a shapefile without any polygons/observations. With ArcGIS Map 10.7 or ArcGIS Pro 2.5, it tells me that I cannot export the shapefile. While ArcGIS aligns the CRS on the fly, it does not allow to compute zonal statistics (because of some mysterious error 99999).
Any suggestions how to solve this problem?


texts: Help with Shortest path Algorithm for enormous bus routes database

texts:
I am building a bus routes map for Chennai. I use Postgis. In our database we have one table having all the routes along with unique ids given for them. Second table comprises of bus stop names with its latitude and longitudes. The third table joins the first and second table based on their ids. We have track records for the routes.

I want to construct a function for finding the shortest route between 2 stops which don’t have a common route. For eg, Two stops A and B exist which don’t have any common routes. Say there are 3 other stops- C,D,E. I want to find out which route should I take to travel the shortest distance. If I leave from A, should I get off at C and from C take a bus to B, or get off at C take another bus to D and from there another bus to B, or go to E and then take a bus. Basically I want to find the shortest path for two stops without a common route. Since my database is enormous that needs to be considered while choosing an Algorithm cause it shouldn’t slow down the process.

Could someone please suggest a suitable algorithm which doesn’t take much time?


texts: How can i find shrines/complex/institute nearby location in qgis and highlight them. The below image is similar to my shape files which contains all information of all attributes,i don't know python, i was using spatial query to find location within area or beside roads, now i want to find complex, malls near to my location and highlight them.


texts: I have two polygons that consist of two separate polygons/lines in QGIS). How can I make this into a multipart polygon, which consists of 1 line in QGIS?


texts:
I am a newbie in Postgres and Postgis and I am trying to do a matching of GPS points to the road network. I need to find closest road from a shapefile and take speed attribute from it. I am using postgreSQL 9.3.7 and postgis 2.1.7.

I have imported a shapefile of the road network into my database. There are many columns, among others: speedlimit, rlid, startdate, enddate, startdistance, enddistance, direction and geom.

I feel unsure about SRID. When I plugged the contents of the .prj file into prj2epsg.org, I got 3006 - SWEREF99_TM as result, and that's what I used to import the shapefile. But I don't know if I need to transform it to 4326, because GPS-points I want to match are longitude and latitude I got from Google maps.

The query:

select distinct(ST_SRID(mytable.geom)) as srid, count(*) from mytable group by srid;

gives:
3006; 2674798

null; 930

I found some solutions and tried them:

Query 1:

SELECT speedlimit, ST_Distance(ST_GeomFromText('POINT(lat long)',3006),geom)
AS distance FROM mytable ORDER BY distance ASC LIMIT 1;

70;6146326.22657711  (this query always gives same result, even if I change GPS points)

Query 2:

SELECT ST_makePOINT(lat long) as gps_point,
ST_Distance( ST_Closestpoint( st_setSRID(r.geom,4326), st_setSRID(ST_makePOINT(lat long),4326) ) ,
ST_makePOINT(lat long) ,true ) as distance_with_c_p, r.speedlimit
FROM mytable r WHERE (ST_Distance_Spheroid (st_setSRID(r.geom,4326),st_setSRID(ST_makePOINT(lat long),4326)
, 'SPHEROID["GRS 1980",6378137,298.257222101]') < 100 )
ORDER BY 2 LIMIT 5

"01010000008849B89047D64C408FFE976BD1222840";1088667.79697218;50
"01010000008849B89047D64C408FFE976BD1222840";1088667.79697218;50
"01010000008849B89047D64C408FFE976BD1222840";3784055.84203355;110
"01010000008849B89047D64C408FFE976BD1222840";5163241.18884849;80
"01010000008849B89047D64C408FFE976BD1222840";5163241.18884849;80

Query 3:

with index_query as (
  select
    st_distance(geom, 'SRID=3006;POINT(lat long)') as distance,
    mytable.speedlimit
  from mytable
  order by geom <#> 'SRID=3006;POINT(lat long)' limit 100
)
select * from index_query order by distance limit 1;

6146326.22657711;70 (this query also gives same result even if I change GPS points)

I don't understand why the nearest distance always is so far away. GPS points are taken from the main roads, and I tested several GPS points.

Is it something about SRID, or what am I doing wrong?


texts:
I have been provided with text (*txt) formatted, tabular climatic data. How can I convert climatic text data to raster in QGIS and ArcGIS? In addition to tabular data, the following information is provided

ncols         480
nrows         600
xllcorner     -20
yllcorner     -59.99999999996
cellsize      0.166667


texts: ArcGIS API for JavaScript - Adding Point(s)

texts: I need to add a single point, I am new to JS and ArcGIS API. I was looking for some code examples on how to create a point and/or a few other points (as needed) I can't use any additional pages/external sources, as in the code needs to be in my single page. I found examples like the button add a point, but I don't need a button I just need to place a point or points with their respective coordinates. Then style as necc. Can someone provide either code or a link so I can do this simple task? Thanks


texts: esri leaflet calcite geojson popup binding

texts:
I am binding a popup to the results of a query, that is working fine. My problem is whenever the map changes (pan, zoom, etc.) the popups no longer work. Is there something I should be doing in a different way?

To reproduce: Select one of the parcel types, click on a popup, pan the map, try again.

http://argis.ualr-iea.org/GFPB/


texts: Arcgis dissolve output has greater area than inputs

texts:
I've been using dissolve to successfully merge polygons (which may or may not overlap) representing mining tenements so that I can calculate their total coverage area (using the geodesic area calculator). However, in one case the output of the dissolve has an 11% greater area than the total area of the input polygons.Have attached a picture of the output.

The input data is from multiple shapefiles. I've selected and merged these into one feature class. I've repaired the geometry of the feature class before running the dissolve.

The dissolve output feature is complex with multiple slivers. I thought that these slivers may be the problem. But eliminating some of them made very little difference to the geodesic area calculation.
I'm baffled!

Does anyone gave any suggestions?


texts: Get the vertices on a LineString either side of a Point

texts:
I have a shapely LineString and defined a shapely Point which lies along the LineString.

How can I find the vertices of the LineString which lie either side of the point? (split the line in two).

This is very similar to this blog post by the author:
http://sgillies.net/blog/1040/shapely-recipes/

Except, instead of cutting by distance along the line, I want to cut given a point that lies along the line.


texts: Adding time band as year to NDVI composite

texts:
I try to add time band to my NDVI time-series and using this function (took it from developers guide), for regression calculation.

It looks like script is working but magnitude of slopes is extremely high.
So, I,d like to add time-band as year, not as a milliseconds, because I am working with yearly composites, how can I to do so?

// This function adds a time band to the image.
var createTimeBand = function(image) {
  // Scale milliseconds by a large constant to avoid very small slopes
  // in the linear regression output.
  return image.addBands(image.metadata('system:time_start').divide(1e18));
};


texts:
I start a new project with base on OpenLayers plugin but if I check the distance on QGIS with distance tool (I choose a straight road) I view 5.466m on QGIS and 4.300m on google (by walking).

How is it possible?

On QGIS 2.6.1 I do not set any CRS


texts: Understanding Google Geocoding API usage restrictions and other free API options?

texts:
I would like to use the google geocoding api to geocode a few thousand addresses.

I understand the daily limit for a free account is 2,500. I want to retrieve the coordinates and then map them in google earth, or maybe in another desktop GIS program.

The usage terms state the geocoding results can only be used in conjunction with a google map.

Does this mean I am out of luck?

I see there are several python geocoding libraries that use the google api, do the same restrictions apply if you are retrieving coordinates using these libraries?

I would like to find a geocoding api that I can use for free to geocode world addresses. From my research, it seems that google has most comprehensive and accurate information, but I want to be sure I can use it.


texts: Copy resolution extent CRS metadata from one raster file to another using raster R

texts:
I have two geotiff raster files, one has all the metadata information and in the other the information was lost. I know that all the metadata information was exactly the same so I want to copy it from one file to the other. I tried to use raster, because I made all the processing in R.

This is the file with metadata

af1_patch <-raster(a_files[6])

af1_patch

class       : RasterLayer
dimensions  : 38400, 38400, 1474560000  (nrow, ncol, ncell)
resolution  : 231.656, 231.656  (x, y)
extent      : -2223901, 6671703, -4447802, 4447802  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs
data source : Forest_patches_AF_1.tif
names       : Forest_patches_AF_1
values      : 0, 255  (min, max)

And this is the file without the metadata

af1_area <-raster(a_files[1])

af1_area

class       : RasterLayer
dimensions  : 38400, 38400, 1474560000  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 38400, 0, 38400  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : africa_AF_1.tif
names       : africa_AF_1
values      : 0, 255  (min, max)

I tried to copy the metadata using:

res(af1_area) <- res(af1_patch)
crs(af1_area) <- crs(af1_patch)
extent(af1_area) <- extent(af1_patch)

but it doesn't work, dimensions and resolution are incorrect and the data values are lost:

af1_area

class       : RasterLayer
dimensions  : 166, 166, 27556  (nrow, ncol, ncell)
resolution  : 53588, 53588  (x, y)
extent      : -2223901, 6671703, -4447802, 4447802  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs

 hist(af1_area)

   Error in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...) :
   cannot make a histogram; need data on disk or in memory


texts: Syntax to view permissions on Oracle table

texts: I want to be able to see what permissions have been granted on a specific object in my Oracle geodatabase. What would my script look like?


texts:
I'm trying to export a 24k image (export set to 1650DPI) using the "export as image" function but I keep getting the "Could not allocate required memory for image" warning.

I don't understand why this is happening as I am running 64bit and have 64gb of RAM installed.

Exporting at a slightly lower DPI of 1400 works and produces a file that is only 50mb. I've checked resource usage and I don't even break the 32gb mark so how can QGIS be running out of memory!?


texts: Elevation profile add-in for ArcGIS not producing graph

texts:
I am trying to create an elevation profile graph using a river line shapefile and the add-in (with ArcMap 10.4.1). So I select my feature, click the "Create a profile for the selected feature" ("Finest" for DEM and blank for "Enter sample distance"). However I get this long error message :

Edit:

File "C:\Users\user\AppData\Local\ESRI\Desktop10.4\AssemblyCache\*\elevationprofile_addin.py", line 152, in GenerateProfile AddGraph(postfix + sid - 1, dem_res, tempWS, tempOutTable, sid, lid, sendToService)
File "C:\Users\user\AppData\Local\ESRI\Desktop10.4\AssemblyCache\*\elevationprofile_addin.py", line 188, in AddGraph arcpy.MakeGraph_management(graph_template, gr, out_graph)
File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 5851, in MakeGraph 'raise e'
RuntimeError: Object: Error in executing tool

It creates two .dbf files and one .cpg file in my scratch folder (also creates a table with 750 lines in my "List by Source" table of contents, but no graph to be found on ArcMap or the scratch folder.

For the record, I also get errors if I try to draw a line or change sampling distance or change DEM res.

Any clues what the problem could be?


texts: Leaflet control search clash with map.on zoomend function

texts:
I would like to know how to solve the clash between Leaflet control plugin and map.on zoomend function

The problem is related to this query:

Displaying selected layers depending on zoom in Leaflet

and as it has been solved, there is still one issue to fix.

I discovered, that once the leaflet search plugin is attached to the map:

  L.control.search({
    layer: L.layerGroup([job2]),
    initial: false,
    propertyName: 'myKey',
    zoom: 18,
    position: 'topleft'
    }).addTo(map);

the map.on('zoomend', function() { doesn't work at the very initial stage (when open or refresh the map to default link).

The problem is shown in the image below:

Under the zoom level 7 all dots are still visible, whereas I set zoomend for zoom 9. When I start zooming in/out on the map, it starts work correctly as perr below:

My code looks as follow:

   L.control.search({
    layer: L.layerGroup([job2]),
    initial: false,
    propertyName: 'myKey',
    zoom: 18,
    position: 'topleft'
    }).addTo(map);

   map.on('zoomend', function() {
   if (map.getZoom() < 9){
if (map.hasLayer(job2)) map.removeLayer(job2);
   }
    else {
     if (document.querySelector("input[name=vm]").checked &&
 !map.hasLayer(job2)) map.addLayer(job2);
   }
  });

When enter the leaflet-search.js plugin, I see the section, which I believe is responsible for this affliction:

   showLocation: function(latlng, title) {  //set location on map from _recordsCache
    var self = this;
    /*
    self._map.once('moveend zoomend', function(e) {

        if(self._markerSearch) {
            self._markerSearch.addTo(self._map).setLatLng(latlng);
        }

    });
    */
    self._moveToLocation(latlng, title, self._map);
    //FIXME autoCollapse option hide self._markerSearch before visualized!!
    if(self.options.autoCollapse)
        self.collapse();

          return self;
      }

Regarding to this I switched off part of the code with: self._map.once('moveend zoomend', function(e) { but the problem remains the same.

Could somebody advice how to remove this issue and make the layers hidden at the initial stage (after refreshing the map)?

Link to my map is here:

https://mariusz-krukar.mkrgeo.pl/en/workshop/media/MDU_demo.html#55.21022,-1.57385,6z

Once open, you will see the dots still visible, whereas zoom level is 6. When yoyu turn into zoom 7 and back to zoom 6 they will be gone until refresh the map.


texts: arcpy; Select by attribute with an If, then statement

texts:
How do I construct an if, then statement based on a select by attribute SQL query or some other method (search cursor?)?

I am asking users to input a unique ID (arcpy.GetParameterAsText(0)), based on that ID I want to search a shapefile field for that value (arcpy.SelectLayerByAttribute_management or arcpy.SearchCursor?), if it is not found, I want to search a different shapefile. When the Unique ID is found I need to export (arcpy.CopyFeatures_management) and isolate the feature.

I am having trouble trying to find something which will indicate if the unique ID is found in the field or not.

I'm essentially looking for something like what occurs with the Verify Expression button event click -> "The expression was verified successfully, BUT NO RECORDS WERE RETURNED".


texts:
I'm trying to merge several layers for a "Count Overlapping Polygon" tool which requires one layer. I have data for fire occurrence polygons in the forest from 1969 to 2016.  I can't get the layers to merge.  Using 10-3-1.  So far, I have tried to simplify a subset of layers in the attribute tables to where there is only 6 columns (FID, Shape, ObjectID, Date, Shape_Leng, Shape_Area). This still doesn't allow them to merged.

The error message is:

  error 001156 : Failed on input OID ; could not write value
  '' to output field .

ArcGIS 10.3 Locate topic
Description

  A value could not be added to the output field. For example, text
  values cannot be added to numeric fields, and text values cannot be
  added to text fields if the values are longer than the field length.

  Solution Change the field type or increase the field length in the
  field map properties.

It suggests to increase the field length but none of the columns have that options.


texts:
Is there any way to represent a polygon geometry that contains a curve accurately in Leaflet (from WKT)? GeoJSON conversion from WKT is insufficient as it digitizes the points of a curve out to an approximation, it's necessary that curves be preserved as their mathematical representation (read:look circular/curved and not like a series of shortly connected lines at lower zoom levels).

I've checked out arc.js and leaflet-omnivore to start with, but my current thinking is that I would have to parse a WKT string into an SVG element and use either the D3 or RaphaelLayer plugins to properly render the layer while preserving curves. I also have performance concerns going the SVG route because I may have many polygons (hundreds to tens of thousands) to render in any given view.

I did find this old question and have looked into the Leaflet.curve and Leaflet.Arc plugins but I don't think either of them could really apply to what I need (I'm also stuck with Leaflet 0.7x for maintainability reasons which eliminates a lot of plugins).

Anyone else have experience with this or a similar problem and have any suggestions or advice on balancing performance/accuracy of the curves?


texts: Python ArcGIS 10.0 Expanding select layer by location from point (walking selection)

texts:
I have an old ArcView 3 tool which I've been tasked with converting to Python. I'm breaking it down into stages, but there's a "walking selection" I can't trap.

To put the query in context, the user selects a site (point data) - selection A. All sites within X kilometres of selection A are then selected (selection B), then all sites within X of kilometres of selection B are selected, and so forth until no new selections are made.

I've tried to use a While... loop with a Get Count to no avail.  The problem I think is because the Get Count will always have a value, so I think I need to do something like Get Count <> Previous Get Count?

#Import ArcPy site package
import arcpy

# Overwrite files
arcpy.env.overwriteOutput = True

 try:
     # Make a feature layer with all the farms
     arcpy.MakeFeatureLayer_management(r"shapefile.shp", "Farms")

    # Make a feature layer containing only the farm of interest
     arcpy.MakeFeatureLayer_management(r"shapefile_select.shp","Farms_select")

    # Apply a selection to farms layer
    fcFarms = arcpy.SelectLayerByLocation_management ("Farms","WITHIN_A_DISTANCE","Farms_select", "2000 Meters")

    # Apply a further selection to farms layer
    fcFarms2 = arcpy.SelectLayerByLocation_management("Farms","WITHIN_A_DISTANCE",fcFarms, "2000 Meters", "ADD_TO_SELECTION")

    while getCount > 0:
        fcFarms3 = arcpy.SelectLayerByLocation_management("Farms","WITHIN_A_DISTANCE", fcFarms2, "2000 Meters", "ADD_TO_SELECTION")
        fcFarms4 = arcpy.SelectLayerByLocation_management("Farms","WITHIN_A_DISTANCE", fcFarms3, "2000 Meters", "NEW_SELECTION")
        getCount = int(arcpy.GetCount_management(fcFarms4))

except:
    print arcpy.GetMessages()


texts: Cesiumjs 1.32 and ArcGIS 10.2 - ArcGisImageServerTerrainProvider is being deprecated - how to provide terrain data

texts:
I tried to connect CesiumJS to ArcGIS by using ArcGisImageServerTerrainProvider.
This works OK if I have a map in simple format BUT:

ArcGIS 10 now supplies data in LERC compression and Cesium can't read it (see http://github.com/Esri/LERC).
From the code I can see that ArcGisImageServerTerrainProvider is deprecated and will be removed in future version.' ArcGisImageServerTerrainProvider` will be removed in 1.32 due to missing TIFF support in web browsers.
see https://cesiumjs.org/Cesium/Build/Documentation/ArcGisImageServerTerrainProvider.html
How can I provide terrain data from ArcGIS 10 to CesiumJS now?


texts: I know there is a setBaseLayer function in OpenLayers 2. Is there a way I can accomplish that in OpenLayers 3?


texts: Splitting address into columns with street name, number and letter

texts:
I have got this column and I want to divide the name, numbers and extra character.

'De waterman 13a' -> 'De waterman' | '13' | 'a'

So I need a column with De waterman a column with 13 and a column with a
and if a field does not have an extra character, it is NULL.


texts: QGIS MSSQL Connection - Raster Data Visible Only as Polygon Footprint

texts:
I've made a connection to ArcSDE on MSSQL from QGIS 3.16.5. I'm able to view the vector datasets in the database and add them to my project view but the raster datasets are only visible as polygon footprints and not the actual raster data.
Does anyone know why I can't view or add raster data from ArcSDE on MSSQL?


texts: comparing raster values over one specific transect

texts:
I have a raster with water temperature data.
I used various kinds of filters over it (smoothing, sharpening filter etc.) in various combination. I want to compare the results, to check which combination is the best for differentiating water spots of different temperatures. I wanted to do cross sections (transects) in few places that interest me and compare differences in values distributions on the line from point A to B.

I looked at various forums and I'm starting to wonder is it even possible?
If anyone knows how to do this (or maybe there is another way)?


texts: ArcPy performance issues with an UpdateCursor when sorting

texts:
I'm hoping that someone will be able to help me here.
I have written a short ArcGIS Python script (mainly cobbled together from other examples) that sorts a feature class ascending based on two fields (route, and sequence). A third field is then populated with an ascending number. An If block restarts the numbering for each route.
The problem is performance. The process is very slow when running on my test data (~8000 records). I plan on running this on 100000 records, so any help or suggestions to improve speed will be appreciated.
Thanks.
James

import arcpy,os
from arcpy import env

inTable = arcpy.GetParameterAsText(0)      #  Input feature class

#  set variables for route and sequence starting values
route_previous = ""
route_current = ""
new_sequence = 1

#  create an update cursor that will append ascending order sequence number to the Service_Type field
#  sorts ascending by Route and Sequence
rows = arcpy.UpdateCursor(inTable,'','','','Route A;Sequence A')

for row in rows:
    route_current = row.Route
#  if the current service location is on the same route as the previous
    if route_current == route_previous:
        row.setValue("Service_Type",new_sequence)
        new_sequence += 1
        rows.updateRow(row)
#  if the current service location is on a different route to the previous
    else:
        new_sequence = 1
        row.setValue("Service_Type",new_sequence)
        new_sequence += 1
        route_previous = row.Route
        rows.updateRow(row)

#  delete the update cursor
del rows


texts: Crop (bounding) box across 180° meridian

texts:
Having a problem in R's raster package -

I can't crop a raster across the 180° meridian - I'm trying to create a bounding box with longitude 120°-230° East

Simple enough, but the following code:

cropbox2 <-c(120,230,70,90)
crop2 <- crop(rt, cropbox2)

Doesn't return anything beyond 180E - it can't seem to cross the 180° meridian

Using minus

cropbox2 <-c(120,-130,70,90)

returns

Error in validityMethod(object) : invalid extent: xmin >= xmax

Surely this is a common problem people have run into?  I'm only trying to crop the N Pacific!